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ABSTRACT 



We characterize the dust in NGC628 and NGC6946, two nearby spiral galax- 
ies in the KINGFISH sample. With data from _3 1 6um to 500/im , dust models are 



strongly constrained. Using the lDraine fc Lil ( 120071 ) dust model, (amorphous sili 



cate and carbonaceous grains), for each pixel in each galaxy we estimate (1) dust 
mass surface density, (2) dust mass fraction contributed by polycyclic aromatic 
hydrocarbons (PAH)s, (3) distribution of starlight intensities heating the dust, 
(4) total infrared (IR) luminosity emitted by the dust, and (5) IR luminosity 
originating in regions with high starlight intensity. We obtain maps for the dust 
properties, which trace the spiral structure of the galaxies. The dust models suc- 
cessfully reproduce the observed global and resolved spectral energy distributions 
(SEDs). The overall dust/H mass ratio is estimated to be 0.0082 ± 0.0017 for 
NGC 628, and 0.0063 ± 0.0009 for NGC 6946, consistent with what is expected 
for galaxies of near-solar metallicity. Our derived dust masses are larger (by up to 
a factor 3) than estimates based on single-temperature modified blackbody fits. 
We show that the SED fits are significantly improved if the starlight intensity 
distribution includes a (single intensity) "delta function" component. We find no 
evidence for significant masses of cold dust (T < 12K). Discrepancies between 
PACS and MIPS photometry in both low and high surface brightness areas result 
in large uncertainties when the modeling is done at PACS resolutions, in which 
case SPIRE, MIPS70 and MIPS160 data cannot be used. We recommend against 
attempting to model dust at the angular resolution of PACS. 

Subject headings: ISM: dust, extinction — ISM: general — galaxies: abundances 
— galaxies: general — galaxies: ISM — infrared: galaxies 



Introduction 



Interstellar dust affects the appearance of galaxies, by attenuating short wavelength 
radiation from stars and ionized gas, a nd contributing IR, submm, mm, and microwave 
emission [for a review, see iDraind ( 120031 )] . Dust is also an important agent in the fluid dy- 
namics, chemistry, heating, cooling, and even ionization balance in some interstellar regions, 
with a major role in the process of star formation. Despite the importance of dust, deter- 
mination of the physical properties of interstellar dust grains has been a challenging task - 
even the overall amount of dust in other galaxies has often been very uncertain. 
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Many previous studies have used far- i nfrare d photometry to estimate the dust properties 
of galaxies. For example, iDraine et al.l ( 120071 ) used global photometry of 65 galaxies in 
the "Spitzer Infrared Nearby Galaxies Survey" (SINGS) galaxies to estimate the total dust 
mass and PAH abundance in each galaxy, and to characterize the intensity of the starlight 
heating the dust. For most of these galaxies the photometry extended only to 160/im, 
although ground-based global photom etry at 850/im was also available for 17 of the 65 
galaxiesn iMunoz-Mateos et al. ( 2009 ) used i mages of the SINGS ga laxies to examine the 
radial distribution of the dust surface density. ISandstrom et al.l ( 120101 1 studied the PAHs in 
the Small Magellanic Cloud SMC on a pixel-by-pixel b asis with a very similar dust model 
as the present work. Very recently, iTotani et al.l (1201 if ) used global photometry in 6 bands 
- 9, 18, 65, 90, 140, and 160/im - to e stimate dust masses for a sample of more t han 1600 
galaxies in the Akari All Sky Survey (llshihara et al.ll2010t lYamamura et al.ll2010l ). In the 
present work, we develop state-of-the-art image processing and dust modeling techniques 
aiming to reliably determine the dust properties in resolved studies. 

The present study makes use of combined imaging by the Spitzer Space Telescope and 
Herschel Space Observatory, covering wavelengths from 3.6/im to 500/im, to produce well- 
resolved maps of the dust in two nearby galaxies. The present study is focused on two well- 
resolved spiral galaxies, NGC 628 and NGC 6946, as examples to illustrate the methodology. 
Subsequent work (Aniano et al. 2 012, in preparation) w ill extend this analysis to all 61 
galaxies in the KINGFISH sample (IKennicutt et al.ll201ll ). 



The paper is organized as follows. In §2] we discuss the data sources. Background 
subtraction and data processing are discussed in ^J3j The dust model is summarized in §4j 
and the fitting procedure is described in £j5j Results for NGC 628 and NGC 6946 are given in 
£j6j The sensitivity of the derived parameters to the set of cameras used to constrain the dust 
models is explored in £J7], where we compare dust mass estimates obtained at high spatial 
resolution (without using MIPS160 or the longest-wavelength SPIRE bands) with estimates 
made at lower spatial resolution (using all the cameras available). We also investigate the 
reliability of the photometry by comparing MIPS70 and MIPS160 images with PACS70 and 
PACS160 images. In $8] we compare dust mass estimates based on spatially-resolved images 
with dust mass estimates based on global photometry, as would apply to distant, unresolved 
galaxies. The "goodness of fit" of different dust models is discussed in §9. The principal 
results are discussed in 510 and summarized in CCD 



Appendices A-D describe the method for image segmentation (i.e., galaxy and back- 



1 Photometry at 850^m was available for 26 galaxies in the SINGS sample, but the data were only reliable 
(and used) for 17 galaxies. 
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ground recognition), background subtraction, estimation of photometric uncertainties in the 
images and estimation of dust modeling uncertainties. Appendix E is a comparison of PACS 
and MIPS photometry. 



Observations and Data Reduction 



NGC 628 and NGC 6946 are part of the SIN GS galaxy sample 



and were imaged by the Spitzer Space Telescope (IWerner et al.ll2004l ). A large subset of the 



Kennicutt et all 120031 ) 



SINGS galax ies are also included in the Herschel Space Observatory Open Time Key Project 
KINGFISH (IKennicutt et al.ll201ll ) and were observed with the Herschel Space Observatory 
jPilbratt et al.lboioh . 



We will use "camera" to identify each optical configuration of the observing instruments, 
i.e., each different channel or filter arrangement of the instruments will be referred to as 
different "camera". With this nomenclature, each "camera" has a characteristic optical 
resolution, spectral response, and point spread function (PSF)@ We will refer to the IRAC, 
MIPS, PACS and SPIRE cameras using their nominal wavelengths in microns: IRAC3.6, 
IRAC4.5, IRAC5.8, IRAC8.0, MIPS24, MIPS70, MIPS160, PACS 70, PACS100, PACS160, 
SPIRE250, SPIRE350, and SPIRE500. For our standard modeling, the PACS and SPIRE 
data were reduced by HIPE followed by Scanamorphos (see Subsection 12 . 2 . 1 1 for details). In 
Appendix E, and Table 6, where we compare PACS data with and without Scanamorphos 
processing, we use PACS (H) 70, PACS (H) 100, PACS (H) 160 to denote PACS data that was 
processed by the HIPE pipeline only. 



Table 1 summarizes the optical resolution of the cameras. 



2.1. Spitzer 



The Infrared Array Camera (IRAC) and the Multiband Imaging Photometer for Spitzer 
(MIPS) cameras on the Spitzer Space Telescope were used to observe all 75 galaxies in 
the SINGS s ample, including NGC 628 and NGC 6946, following the observing strategy 
described by IKennicutt et al.l (120031 ). Spectroscopic observations of selected regions were 
also obtained, although not used in the current study. Aniano, Reyes, Draine et al (2012, 



2 For example, the 70/imand 100/im channels of the PACS instrument use the same optical path in the 
telescope, but differ in the filter used, and thus have different PSFs, so wc will consider them as different 
"cameras" . 
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in preparation; hereafter ARD12) use the spectroscopic observations of selected regions to 
further constrain the dust modeling. 



2.1.1. IRAC 



IRAC ( jFazio et al.l 120041 ) imaged the galaxies in 4 bands, centered at 3.6/im, 4.5/xm, 
5.8/^mand 8.0/im. The images were processed by the SINGS Fifth Data Delivery pipelined 
The IRAC images are calibrated for point sources. Photometry of extended sources requires 
so-called "aperture corrections" . We multiply the intensities in each pixel by the asymptotic 
(infinite radii) value of the aperture correction (i.e., the aperture correction corresponding 
to an infinite radius aperture). We use the factors 0.91, 0.94, 0.66 and 0.74 for the 3.6/xm, 
4.5/xm, 5.8/xmand 8.0/im bands, respectively, as described in the IRAC Instrument Hand- 
book (V2.0.lfl 



2.1.2. MIPS 



described in Kennicutt et al 



Imaging with MIPS flRieke et al.ll2004h was carried out following the observing strategy 



(ESS). The data were reduced using the LVL (Local Volume 



Legacy) project pipeline^ A correction f o r non linea rities in the MIPS70 camera was applied 



by the team, as described by iDale et al.l (120091 ) and iGordon et al.l (120111 ). 



2.2. Herschel 

The Photodetector Array Camera and Spectrometer for Herschel (PACS) and the Spec- 
tral and Photometric Imaging Receiver (SPIRE) cameras on the Herschel Space Observatory 
are being used to observe the 61 galaxies in the KINGFISH s ample, in particu l ar NG C 628 



and NGC 6946, following the observing strategy described by lKennicutt et al.l (120111 ). The 
maps were designed to cover a region out to > 1.5 times the optical radius R op t, with good 
signal-to-noise (S/N) and redundancy. The depth of the PACS images at 70/xmand 160/imis 
less than that of MIPS, but the higher resolution of PACS is able to better single out compact 



3 Details can be found in the data release documentation, http://data.spitzer.caltech.edu/popular/sings/20070410-cnhance< 



http://irsa.ipac.caltcch.edu/data/SPITZER/docs/irac/iracinstrumcnthandbook/IRACJnstrumcnt_Handbook.pdf 



5 Dctails can be found in the data release documentation, http://data.spitzer.caltech.edu/popular/lvl/20090227-enhanced/( 
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star-forming regions. 



2.2.1. PACS 



NGC 628 and NGC 6946 were observed with the PACS instrument on Herschel (IPoglitsch et al. 



201(f ) on 2010 Jan. 28 (NGC 628) and Mar. 10 (NGC 6946), using the "Scan Map" observa- 
tion mode. The PACS images were first reduced to "level 1" (flux-calibrated brightness time 
series, with attached sky coordinates) using HIPE (jOttll2010l ) version 5.0.0, and ma ps ( "level 
2") were created using the Scanamorphos data reduction pipeline (IRoussell 120121). version 



16.9. This reduction strategy includes the latest PACS calibration available ( IMiiller et al. 
201ll ). and aims to preserve the low surface brightness diffuse emission. 



Additionally, PACS data was reduced completely (i.e., to "level 2") using HIPE. We 
used HIPE version 5.0.0, and further divide the fluxes by 1.119, 1.151, and 1.174 for the 
cameras PACS70, PACS100, and PACS160 respectively to account for the latest calibration 
of the cameras. Both PACS data reduction strategies present strong discrepancies with the 
corresponding MIPS photometry, as is discussed in §|F] The HIPE pipeline removes a large 
fraction of the flux in the low surface brightness areas, and in our work it was only used for 
comparison, i.e., we found (Appendix F) that the Scanamorphos pipeline is more reliable, 
and we only employ the Scanamorphos reductions when PACS data are used in the dust 
modeling. 



2.2.2. SPIRE 



The 2 galaxies were observed with the SPIRE instrument (iGriffin et al.ll2010l ) on 2009 
Dec. 31 (NGC 6946) and 2010 Jan. 18 (NGC 628). The data were first reduced to "level 1" 
using HIPE version spire-8. 0.3287, followed by Scanamorphos version 17. C0. The assumed 
beam sizes are 435.7, 773.5, and 1634.6 arcsec 2 for SPIRE250, SPRE350, and SPIRE500, re- 
spectively. Additionally, we excluded discrepant bolometers from the map, and adjusted the 
pointi ng to match the MIPS24 map. Data reduction details can be found in lKennicutt et al. 

(honh . 



6 Sec also 



ftp://hsa.esac.esa.int/URDj-ep/KINGFISH.DRl_SPIRE/ 
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2.3. Atomic and molecular gas 



NGC 628 and NGC 6946 are in The H I Nearby Galaxy Survey (THINGS; IWalter et al. 
20081 ). with ~6" resolution H I 21-cm imaging by the NRAO Very Large Array. The H I 



21-cm emission is assumed to be optically-thin, in which case the H I surface density is 
directly proportional to the 21-cm line intensity. 

The CO J = 2 — > 1 transition was observ ed with ~13" ang ular resolution by the HERA 



CO Line Extragalactic Survey (HERACLES; iLeroy et al.l 120091 ) using the HERA multipixel 
receiver on the IRAM 30-m telescope, with estimated uncertainties of ±20%. 

As is usual, the H 2 mass surface density is taken to be proportional to the CO J = 2 — > 1 
line intensity. Xco is the ratio of H 2 column density to CO J = 1 — > intensity integrated 
over the line profile. In what follows, we define 

Xco / -i \ 



X 



and w e assume a J = 2 — > 1/ J 



20091 ). 



= 10 20 H 2 cm- 2 (Kkm s" 1 )" 1 ' 
1 — > antenna temperature ratio i? 2 i 



0.8 ( Leroy et al. 



The conversion factor Xco is uncertain. Xcq [2 q ~ 2 is the value normally adopted for 

molecu lar gas in the Milky Way (IDame et al.ll200ll ; lOkumura et al.ll2009f). iPlanck Collaboration et al 
fepllbh found X co ,2o = 2.54 ± 0.13 for the Milky Way, braine et all toor\) found that 
Xco, 20 ~ 4 appeared to g ive the most reasonable dust/H mass ratios for the SINGS galaxy 
sample. iBlitz et al.l (120071) fou nd Xco, 20 ~ 4 to be the best overall value for galaxies in the Lo- 
cal Group. iLerov et al.l (boilh found X co ,2o = 1.2 - 4.2 for M 31, M 33, and the Large Mag- 
ellanic Cloud (LMC), and the very high values X co ,2o = 14 and X co ,2o = 32 for NGC 6822 
and the Small Magellanic Cloud (SMC). In the present work, we take Xco, 20 = 4 ± 1 to be 
the best overall value for NGC 628 and NGC 6946. In §6] (Fig. |3]) we show the H gas maps 
and dust/H mass ratio for Xco,2o = 2, 3, 4. 

The H I and H 2 masses are added to generate maps of the total H surface density Eh- 
Including Helium would give S gas « 1.38 x E H , but in our discussion will use S H as it is the 
"observable" , avoiding uncertainties in the Helium abundance. 



3. Image Analysis 

Before modeling the spatially resolved SEDs, it is necessary to adjust the images in 
a number of ways to ensure meaningful results. In the following sections we describe the 
image analysis steps in detail. These steps include background estimation and subtraction, 
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convolution to a common PSF, and resampling of the convolved images to a common pixel 
grid. We also correct the final images for missing or bad data in the original images and 
estimate the uncertainties on the flux in each pixel. The pixel sizes used in the final maps 
are chosen to Nyquist-sample the final-map PSFs. 

All camera images are first rotated so that North is up, and then trimmed to a common 
sky region. We then estimate and remove a background from each image, as described in 
Appendix [B] Following this, we convolve the images to a common point spread function 
(PSF), and resample all the images on a common final- map grid. We correct the final 
images for the bad pixels (or missing data) in the original images, as described in Appendix 
ICl Finally we use the background pixel dispersion to estimate the pixel flux uncertainties, 
as described in Appendix |D] 

3.1. Background recognition and Subtraction 

We use a multistep algorithm to generate a Background Mask, consisting of the area not 
covered by either the target galaxy or other discrete sources (e.g., recognizable background 
galaxies or foreground stars). This algorithm is described in detail in Appendix |A] and avoids 
overestimating the background in cameras with low S/N. 

After the background mask is generated, for each image we estimate the best-fit back- 
ground "tilted plane", as described in Appendix [A] After the final best-fit background "tilted 
plane" has been found for each camera, it is subtracted from the original images. The back- 
ground subtraction is performed on each original map independently, using its native pixel 
grid. 

The result is, for each camera, a background-subtracted image with its native pixel grid 
and PSF. Figures [1] and [2] show the resulting background-subtracted images. 

3.2. Convolution to a common PSF 

In order to perform any resolved dust study, it is necessary to convolve all the images 
to a common PSF. To generate maps with appropriate wavelength coverage to perform the 
dust modeling, the natural final-map PSFs to use are those of the PACS160, SPIRE250, 
MIPS70, SPIRE350, SPIRE500, and MIPS160 cameras. For a given final-map PSF, only a 
subset of cameras may be transformed into it reliably, and we proceed to investigate the most 
reasonable compatible camera combinations, considering the tradeoff between (1) angular 
resolution and (2) availability of long-wavelength data to constrain the dust models. After 



- 9 - 



0.6 - Sorry, this plot is 



04 too large for this preview.. 



Sorry, this plot is 



too large for this preview.. 



0.6 - Sorry, this plot is 



04 too large for this preview.. 



Sorry, this plot is 
too large for this preview. 



Sorry, this plot is 



too large for this preview.. 



Sorry, this plot is 



too large for this preview. 



Sorry, this plot is 



too large for this preview.. 




Sorry, this plot is 



too large for this preview.. 




24.2 24.1 
RA (deg) 



24.2 24.1 
RA (dcg) 



Fig. 1. — Background-removed images of NGC 628. Top row: IRAC3.6, 5.8, 8.0 (Left, middle, and right 
row respectively). Second row: MIPS24, 70, 160. Third row: PACS70, 100, 160. Bottom row: SPIRE250, 
350, 500. The white contour is the boundary of the galaxy mask, within which the data allow reliable 
estimation of the dust properties. 
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Fig. 2.- 



Same as Figure [U but for NGC 6946. 
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choosing the appropriate PSF for a given set of cameras, we transform all th e background- 



subtra cted images to this common PSF using convolution kernels described by lAniano et al. 



( 120111 ). In §6 we will focus on maps generated at three final-map PSFs: PACS160, SPIRE250, 
and MIPS 160. PACS160 is the PSF with smallest FWHM (full width at half maximum) 
that allows use of enough cameras to constrain the dust SED (IRAC, MIPS24, and PACS70, 
100, 160). SPIRE250 allows use of the same cameras as PACS160 plus the SPIRE250 
camera. Adding the 250/xm constraint produces more reliable maps. The MIPS160 PSF 
allows inclusion of all the cameras (IRAC, MIPS, PACS, SPIRE), therefore producing the 
most reliable dust maps; this will be our "gold standard" . 

The convolution kernels assume that the PSFs can be approximated by rotationally 
symmetric functions. In general, a convolution kernel will relocate flux in the images to 



transform them to a desired PSF. lAniano et al.l (120111 ) developed a criterion for camera 



compatibility, to determine which cameras can be reliably transformed into a given PSF. 
Essentially, a PSF can be safely transformed into another PSF with similar extended wings 
provided that the final FWHM is larger than the original. When the extended wings of 
the two PSFs are dissimilar, the criterion involves quantifying the amount of energy that 
a kernel should remove from the extended wings of the first PSF, as this power removal is 
correlated with the risk of introducing artifacts in the convolved image. The performance of 
the convolution kernels is excellent: for "safe" pairs of PSFs, the discrepancies between the 
convolved narrower PS F and the broader P SF is smaller than the uncertainties in determining 
the PSFs themselvefl U m h «t nil implemented a method to construct convolution 



kernels for non rotationally symmetric PSFs, using the theoretical PSFs for the Herschel 
cameras. This method relies on the theoretical PSFs, whereas in the present method we are 
able to use rotational averaging to empirically characterize the extended wings of the actual 
PSFs measured using observations of saturated point sources. Table [TJ lists the resolutions 
of the cameras, the pixel size in the final-map grids used, and the other cameras that can be 
used at this resolution. 

The CO J = 2 — > 1 maps (used to generate the H 2 maps) are provided in rotation- 
ally symmetric gaussian PSFs with 13.4"FWHM, which can be safely transformed into 
all the final-map PSFs used. The original H I maps have non-circular gaussian PSFs. 
(FWHM=9.30" x 11.88" for NGC 628, and FWHM=5.61" x 6.04" for NGC 6946). When 
convolving the H I maps into the final map resolutions, we will use kernels generated for 
rotationally-symmetric gaussian PSFs with FWHM=^/9.30 x 11.88 " and \/5.61 x 6.04" for 
NGC 628 and NGC 6946, respectively. 



7 Table 4 of lAniano et al 
have D < 0.064. 



(|201ll) quantifies the corresponding PSF mismatch: all the kernels employed 



- 12 - 



Table 1: Image resolutions 





FWHM a 


50% power" 


Final grid 


Compatible 


Camera 


(") 


diameter (") 


pixel 6 (") 


cameras c 


IRAC3.6 


1.90 


2.38 




not used as a final-map PSF 


IRAC4.5 


1.81 


2.48 




not used as a final-map PSF 


IRAC5.8 


2.11 


3.94 




not used as a final-map PSF 


IRAC8.0 


2.82 


4.42 




not used as a final-map PSF 


PACS70 


5.67 


8.46 




not used as a final-map PSF 


MIPS24 


6.43 


9.86 




not used as a final-map PSF 


PACS100 


7.04 


9.74 




not used as a final-map PSF 


PACS160 


11.2 


15.3 


5.0 


IRAC; MIPS24; PACS 


SPIRE250 


18.2 


20.4 


8.0 


IRAC; MIPS24; PACS;SPIRE250 


MIPS70 


18.7 


28.8 


10.0 


IRAC; MIPS24,70; PACS; SPIRE250 


SPIRE350 


24.9 


26.8 


10.0 


IRAC; MIPS24,70; PACS; SPIRE250,350 


SPIRE500 


36.1 


39.0 


15.0 


IRAC; MIPS24,70; PACS; SPIRE 


MIPS160 


38.8 


58.0 


16.0 


IRAC; MIPS; PACS; SPIRE 



a Values from Aniano et al. (2011 ) for the circularized PSFs. 

b The pixel size in the final-map grids is chosen to Nyquist-sample the PSFs. 

c Other cameras that can be convolved into the camera PSF (see text for details). 
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3.3. Uncertainty estimation 

For each camera, after the image processing (rotation to RA-Dec, background subtrac- 
tion, convolution to a common PSF, and resampling to the final grid) the flux in each final 
pixel is a (known) linear combination of the flux of (in principle) all the original pixels of 
the camera. If the statistical properties of the uncertainties in the original pixel fluxes were 
known, it would be possible to propagate these uncertainties (and their statistical proper- 
ties) to each final pixel. The original maps oversample the beam, and have artifacts that 
extend over several pixels, so realistic statistical properties of the uncertainties are difficult 
to determine. We therefore estimate the uncertainties directly in the final (post-processed) 
image. 

Using the pixels of the background mask (adapted to the final- map grid), we measure the 
dispersion of the background pixels (which includes noise coming from unresolved undetected 
background sources, image artifacts and detector noise) as described in Appendix [D] By 
comparing the MIPS and PACS images, we can also estimate a calibration uncertainty, as 
described in Appendix |D] 



4. Dust Model 

The composition of interstellar dust remains uncertain, but models based on a mixture 
of amorphous silicate grains and carbonaceous grains have proven successful in reproducing 



the m ain observed properties of interstellar dust. We employ the dus t model of iDraine fc Li 



( 20071 . hereafter DL07), using "Milky Way" grain size distributions ( IWeingartner &: Draine 



200ll ). The DL07 dust model has a mixture of amorphous silicate grains and carbonaceous 



grains, with a distribution of grain sizes, chosen to rep roduce the wavelength depe ndence 



of interstellar extinction within a few kpc of the Sun ( IWeingartner &: Draind l200ll ). The 
silicate and carbonaceous content of the dust grains was constrained by observations of the 
gas phase depletions in the interstellar medium. 

The bulk of the dust in the diffuse ISM is heated by a general diffuse radiation field 
contributed by many stars. However, some dust grains will happen to be located in regions 
close to luminous stars, such as photodissociation regions (PDRs) near OB stars, where the 
starlight heating the dust will be much more intense than the diffuse starlight illuminating 
the bulk of the grains. Since our pixels have a large physical size (~500 pc side for MIPS 160 
PSF), we will assume that, in each pixel, there is dust exposed to a distribution of starlight 
intensities. 
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4.1. Carbonaceous Grains 

The carbonaceous grains are assumed to have the properties of polycyclic aromatic 
hydrocarbon (PAH) molecules or clusters when the number of carbon atoms per grain Nq < 
10 5 , but to have the properties of graphite when Nq 3> 10 5 , with an ad-hoc smooth transition 
between the two regimes. 

A carbonaceous particle of equivalent radius a is taken to have absorption cross section 

C a bs(a, A) = (1 - / g )A r c o r PAH(H : C, A) + / g C ab s (graphite, a, A) , (2) 

where <Tpah(H:C, A) is the absorption cross section per C for PAH material with given H:C 
ratio, C a bs (graphite, a, A) is the absorption cross section calculated for randomly-oriented 
graphite spheres of radius a, and the graphite "weight" is taken to be 

/ g = 0.01 for a<a c (3) 
= 0.01 + 0.99 [1 - [ajaf] for a > a c , (4) 

where we take the transition radius a c = 0.0050/im. Carbonaceous grains with a > a c are, 
therefore, treated as having a graphitic component, with the graphite weight f g — > 1 for 
a 3> a c . However, with f g = 0.01 for a < a c , even very small PAHs are assumed to have a 
small continuum opacity underlying the PAH features. 

The PAH absorption cross section can be represented as the sum of a num ber of vibra- 



tional features with specified central wavelength, FWHM, and band strength. iLi fc Draine 



(120011 . hereafter LD01) presented a s et of resonance pa rameters that appeared to be con 



sistent with pre-Spitzer observations. ISmith et all (120071 ) used the PAHFIT fitting software 



with the SINGS spectra to improve observational determinations of central wavelengths, 
shapes, and overall strengths of the PAH emission profiles; DL07 used these results to adjust 
the LD01 profile parameters. Subsequent modeling of the SINGS nuclear spectra (ARD12) 
led to some additional small changes in some of the PAH band strengths. In the present 
study we employ the PAH cross sections from DL07 and ARD12. 



Draine fc Lil (120071 ) adopted the model put forward by iDraine fc Led (119841 . hereafter 
DL84) for the far-infrared properties of graphite. Graphite is a highly anisotropic material, 
with very different responses for E || c and £lc, where the c axis is normal to the "basal 
plane". DL84 included "free-electron" contributions 5e^_, 5e^ to the dielectric tensor, using 
a simple Drude model for the free electron response, 



(5) 
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where oj = 2-kc/X, and to allow for size effects the "mean free time" r is taken to be 

r_1 = r buik + u F/a , (6) 

where Tb u ik is the mean free time in bulk material, vf is the Fermi speed, and a is the grain 
radius. 

For i? 1 ewe continue to use the graphite dielectric function from DL84. However, for 
E || c, the DL84 free electron model for Sej resulted in an opacity at A > 100/xm that gav e 
somewhat more emission than observed in the Milky Way cirrus by iFinkbeiner et al.l ( 119991 ) . 
In addition, the free-electron model used by DL84 for E || c produced an opacity peak 
near 33/xm that does not give a good match to InfraRed Spectrograph (IRS) observations of 
emission from regions where the grains are hot enough to radiate near 33/xm. To broaden 
the opacity peak, we now take the free-electron contribution for E || c to be 



5el\ 



OJ] 



urn 



+ iUTj 



(7) 



with oj. 



and oj. 



),2 



The Tj are obtained from eq. ([6]) with 
0.88 x 10~ 14 s. This gives a d.c. electrical conductivity 
j{^p jTj)/4ir = 5.62 x 10 13 s _1 = 62.5 mh o cm" 1 , with in the range reported for high-quali ty 



1 x 10 14 s _1 , 

-14 , 



T bu ik,i = 3.51 x 10" "s, r bulk , 2 
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2 x 10 14 s _1 . 
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graphite crystals at 300 K [~ 1 mho cm -1 (jKleinlll962l ) to ~ 200 mho cm -1 f lPrimakill956[ )]. 
This d.c. conductivity is larger than the value 30 mho cm -1 adopted by DL84; the increased 
conductivity lowers the FIR emission and br ings the overall emission spectrum into better 
agreement with the observed spectrum from IFinkbeiner et al.l (119991 ) . We take vp = 3.7 x 
10 6 (1 + T/255 K) 1 ^ 2 cms -1 . The two-component free-electron form of eq. ([7]) is not intended 
to have physical significance. It is adopted because it is analytic, satisfies the Kramers- 
Kronig relations, gives a reasonable value for the d.c. electrical conductivity, and results in 
an opacity that is less peaked than the original single- component form §5$). The resulting 
graphite opacity varies as A~ 2 for A > 200/im. 



4.2. PAH Abundance ^pah 

As discussed above, the PAHs are part of the carbonaceous grain population. The 
PAH abundance is measured by the parameter ^pah; defined to be the fraction of the total 
grain mass contributed by PAHs containing iVc < 10 3 C atoms. The PAH size distribution 
used in the DL07 models extends up to PAH particles containing iVc > 10 5 C atoms (a > 
6.0 x 10 _7 cm). However, IRAC photometry at 5.8yumand 8.0/imis sensitive primarily to 
PAHs with A^c ^ 10 3 C atoms, small enough so that single-photon heating can result in 
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significant 8/im emission (see, e.g., Fig. 7 of iDraine Sz Lil 120071 ). For the size distribution 



in the DL07 models, the mass fraction contributed by PAH particles with Nq < 10 6 is 
1.478 gpAH- 

WD01 constructed grain size distributions with different values of (/pah that were com- 
patible with the average extinction curve in local diffuse clouds. Such models are possible 
for gpAH ^ 0.046, with part of the 2175A extinction feature contributed by the PAHs, and 
part contributed by small graphitic grains. For (/pah ^ 0.046 the predicted 2175A feature 
from the PAHs alone would be stronger than the average observed 2175A feature in local 
diffuse clouds 



Nevertheless, because the PAH abundance in other regions could conceivably exceed 
the value in the local Milky Way, the WD01 dust models have been extended by simply 
adding PAHs to the (/pah = 0.046 model, with no adjustment to the populations of silicate 
or larger carbonaceous grains. These models produce stronger emission in the PAH emission 
features, particularly in the IRAC8.0 band. The models were extended to (/pah = 0.10. The 
models with (/pah > 0.046 have a 2175A feature strength larger than the average value in 
the local ISMjj The model set was also extended down to (/pah = (the smallest value 
of (/pah considered by WD01 was 0.0047). Models were computed in a grid of (/pah = 
0, 0.01, 0.02, ... 0.10, and linearly interpolated to a grid with spacing A(/pah = 0.001. 



4.3. Amorphous Silicate Grains 



In the DL84 dust model, the amorphous silicate absorption in the infrared was modeled 
by a set of damped Lorentz oscillators, resulting in an opacity varying as A~ 2 for A 3> 
25/im. However, the COBE-F IRAS measurements of the A > llO^m emission spectrum o f 
dust at high galactic latitudes (j Wright et al. 1991 ; Reach et al. 1995 ; Finkbeiner et al. 1999 ) 
were not ac c uratel y reproduced by the A -2 opacity of the DL84 graphite-silicate model. 
Li fc Draind (|200l[ ) therefore made an ad-hoc modification to Im(e) for amorphous silicate 
at A > 250/mi, so that it is no longer a simple power-law [with the corresponding changes to 
Re(e) required by the Kramers-Kronig relations] . The adj ustments to Im(e) were not large - 
the modified amorphous silicate Im(e) (ILi &: Draind l200ll ) is within ±12% of that for DL84 



iDraine et al.l (|2007l ) in a global study of 61 galaxies in the SINGS sample found a median value of 
gPAH = 0.034. 

In our modeling, we find best-fit qpah < 0.05 for virtually all of the galaxy pixels, compatible with the 
average observed 2175A feature in local diffuse clouds, i.e., models with qpah > 0.05 are not needed, except 
for error estimation. 
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amorphous silicate for A < llOOyum - but these modest adjustments brought the emission 
spectrum for the dust model into fairly good agreement with o bservations of the emission 
spectrum of high-latitude dust (see Fig. 9 of iLi fc Draindl200ll ). The adopted opacity has 
no dependence on the grain temperature T. 

It is possible that the amorphous silicate opacity may in actuality be T-dependent 
( Menv et all 120071) . and s ome authors have argued that this is indicated by ob servations 
JParadis et al.lboiol . [201lh . The "two-l evel-svstem" m odel of lMeny_el_alJ J2007L with the 
standard parameters recommended by iParadis et all (1201 ll ). has the far-infrared spectral 
index /3 = din k/ din v near A = 500yum varying from ~ 2 to ~1.3 as T increases from 
10K to 50K. However, in the present study we find that the DL07 dust model is able to 
satisfactorily reproduce the observed spatially-resolved SEDs, as well as the global emission. 
At least for near-solar metallicity galaxies such as NGC 628 and NGC 6946, dust models 
with T-dependent opacities do not appear to be required. 



4.4. Dust Heating 

Each dust grain is assumed to be heated by radiation with energy density per unit 
frequency 

u „ = Ux u^ MP83 (8) 
where U is a dimensionless s caling factor and 

M MMP83 ig 

the interstellar radiation field (ISRF) 



estimated by iMathis et al.1 (Il983[ ) for the solar neighborhood. We ignore variations in the 
spectral shape. 

Each pixel in our modelling will be larger than 2 x 10 4 pc 2 , so it will contain ISM in a 
variety of physical enviroments. A fraction (1 — 7) of the dust mass is assumed to be heated 
by starlight with a single intensity U = £/ mm (i.e., heated by a diffuse ISM radiation field), 
while the remaining fraction 7 of the dust mass is exposed to a power-law distribution of 
starlight intensities between U min and U max with dM/dU oc U~ a . 

The starlight heating intensities are thus characterized by four parameters: 7, Umin, 
U max , and a, where the fractional dust mass dM A {U) heated by starlight intensities in (U, U + 
dU) is 

(1 - i)S(U - U min ) + T -(g— il— U~ a for U mm <U< U m (9) 



M d ,tot \dUj- - 
for a 7^ l,and 



1 fdM, 



id 



(1 - <y)6{U - U min ) + 1 -— -—-U~ l for U mm <U< U max (10) 



M d>to t \ dU J In (U max /U min ) 
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for a = 1, where Mdtot = f^ max (dM&l dU')dU' . More complicated starlight heating distri- 

J ^min 

butions could be contemplated, but we find that the simple 4-parameter (7, t/ min , C/ max , a) 
model of eq. (P - DUD appears able to usually provide an acceptable fit to observed SEDs in 
star-forming galaxies with near-solar metallicities. 

Galliano et al (2011) recently claimed that the emission from the Large Magellanic 
Cloud (LMC) can be reproduced using a starlight distribution function that lacks the "delta 
function" component of eq (M - [TOT) . i.e., fixed 7 = 1. However, we show in §21 that adding 
the "delta function" component significantly improves the quality of the fit for the pixels in 
NGC 628 and NGC 6946. 

Many authors choose to fit the A > 70 /im emission using a blackbody -B„(Td) multiplied 
by a power law opacity oc A -/3 . The best-fit value of is closely related to our heating 
parameters U min , 7, a. In a subsequent work (Aniano & Draine 2012, in prep.) we show 
that Td ~ 20 U ^h? K, when the DL07 SED is approximated by a blackbody multiplied by a 
power law opacity. 

Given that there may be significant regional variations in the starlight spectrum, U 
should be interpreted not as a measure of the starlight energy density, but rather as the 
ratio of the actual dust heating rate to the heating rate for the MMP83 radiation field. 

The fraction of dust luminosity emer ging in the PA H features does depend on the 



spectrum of the starlight heating the dust. iDraind fj2011al ) showed that the fraction of the 



dust emission appearing at 8/im increases by a factor of 1.57 as the starlight spectrum is 
changed from MMP83 to a 20 kK blackbody (cut off at 13.6eV). If the actual hv < 13.6eV 
starlight spectrum is harder (softer) than the MMP83 spectrum assumed in the models, we 
will overestimate (underestimate) ^pah- 



4.5. Contribution of Direct Starlight 



Starlight enters in the dust modeling in two ways: via dust heating (as discussed in 
§4.4p . and as a direct starlight component (i.e., direct starlight escaping the observed region). 
Our main goal in the present work is to study the properties of the dust and the starlight 
heating the dust, s o we adopt a simple mod el for the direct starlight component. Following 
Bendo et al.l (120061 ) and iDraine et al.l (120071 ). we approximate the A > 3/xm stellar emission 
from the galaxy as simply 

f*(a) = n*B„(A,r*) , (11) 

where is the solid angle subtended by the stars, B v is the blackbody function, and 
T* = 5000 K is a representative photospheric temperature to approximate the integrated 
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stellar emission at A > 3yum. The direct starlight contribution is thus adjusted by only the 
parameter Q+. Direct starlight will only contribute significantly to the IRAC bands, and very 
marginally to MIPS24. In ARD12 corrections arising from photospheric absorption as well 
as emission from hot circumstellar dust around AGB stars are studied. These corrections 
are only important in elliptical galaxies with little interstellar medium, and the results in the 
present paper would be virtually unchanged if included. We also neglect possible reddening 
at the wavelengths (A > 3.6/im) in the present study. 

Although in principle the direct starlight and heating starlight parameters should be 
connected, the uncertain and complex distribution of dust and stars within the galaxies make 
such connection very complex^!. Our simplified treatment of the direct starlight should only 
be regarded as a way of "removing" the direct starlight component from the near infrared 
photometry so we can have an estimate of the dust emission. 



4.6. Dust Model Emission 

For each given set of dust parameters (gpAH, 7, ^min, ^max, «), and the given chemical 
composition, grain size distribution and grain properties, the dust emission spectrum is 
computed from first principles. 

First, for each given starlight heating parameter U, the temperature distri bution of the 



dust grains (including the PAH component) is computed as described elsewhere (IDraine fc Li 



200ll . 120071 ). This is performed for a logarithmically-spaced grid of 41 U values from 0.01 to 
10 8 . From the temperature distribution functions, model spectra are computed and stored. 
To obtain spectra for intermediate U values, we interpolate. 

Secondly, for each starlight heating distribution, the specific power spectrum per unit 
dust mass p u (model) is computed. We essentially have two independent heating starlight 
intensity distributions, the "delta function component" (i.e., the dust exposed to U — £/ min ) 
and the "power law component" (i.e., the dust heated by starlight with U min < U < U max ). 
Lastly each p v (model) is convolved with the various spectral response function^ to ob- 
tain the predicted photometry (p(model))fc for each camera k, with nominal wavelength 
Afc. We construct a library of model emission for a finely-sampled grid of parameters 



10 The mean stellar brightness emerging from each pixel is not the same as the starlight intensity seen by 
the dust grains, and therefore we do not expect a tight correlation between the dust heating and f2+ 

11 We employ the appropriate spectral response function for each camera. For SPIRE, we use the spectral 
response functions appropriate for extended sources, as described in the SPIRE observers manual. 
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'<?pah, U min , U max , a) and (<?pah, U n 



5. Determining the Dust and Starlight Heating Parameters. 

NGC 628 and NGC 6946 are well resolved, with each galaxy providing many independent 
pixels, even at MIPS160 resolution. For each pixel j, we find the model of dust and starlight 
that best reproduces the observed SED, within the modeling scheme described by DL07. 

As discussed in §4.5[ starlight enters the fitting in two ways: via direct starlight in the 
pixel, and by heating the dust. The direct starlight contribution to pixel j is adjusted by 
varying only the parameter The heating starlight intensity is characterized by four pa- 
rameters: 7j, Urainj, U m3iXt j, and aij, where the dust mass dM& heated by starlight intensities 
in (U, U + dU) is given by equation d^J). For each pixel j we adjust the total dust mass Mdj, 
the PAH abundance parameter qpAKj (PAH mass fraction), and the characteristics of the 
starlight heating the dust in that pixel. If mid-IR photometry is unavailable, one loses the 
ability to constrain gpAH, but (adopting some arbitrary value of (/pah for the modeling), the 
dust mass estimation itself would be largely unaffected. Fortunately, we do have mid-IR 
coverage of our galaxies, so we can obtain full qpah maps for them. The present modeling 
assumes the grains to be heated by a standard starlight spectrum (corresponding to the 
starlight in the local ISM), and to have a fixed balance between PAH neutrals and ions. 
When fitting to IRAC, MIPS, PACS, and SPIRE photometry, the best-fit value of q PA u is 
then essentially proportional to the strength of the (nonstellar) IRAC8.0 band power relative 
to the total IR power. 

We will find that jj <C 1 in nearly all regions where the dust luminosity surface density 
Sl c1 j > 10 7 L Q kpc~ 2 : here, U m i n j is presumed to represent the diffuse interstellar med ium, 



or the counterpart to the "infrared cirrus" component of the galaxy (lLow et al.lll984f ). ac- 
counting for the bulk of the dust mass in pixel j. The small fraction jj of the dust mass 
exposed to starlight intensities U > U inin j is presumed to correspond primarily to dust in 
star-forming regions. 



The model flux density in camera k is 

4vrD 2 

where {F*(\))k is the direct contribution of starlight given by Eq. ( TTTT) convolved with the in- 
strumental response function. The dust model is characterized by M d , (Jpah, 7, U min , ?7 max , «}• 



F(model, A fc ) = (F*(A)) fe + 7 -^(p(model)) fc , (12) 



In principle, t/ maxJ could be treated as an adjustable parameter. Previous work (IDraine et al. 



20071 ) has shown that the quality of the global fit to the SED is relatively insensitive to the 
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choice of U max . We experimented by allowing U max j to be fitted I in the resolved maps, and 
found that the best value for most of the pixels is t/ maxJ = 10 7 . In the pixels where the 
resulting best-fit value is not 10 7 , fixing £/ max j = 10 7 did not decrease the quality of the fit 
significantly, i.e. the total \ 2 is essentially the same, as shown in §9. Allowing U max to be 
fitted in the range 10 3 < U max < 10 7 or fixing it to U m3X = 10 7 does not produce appreciable 
changes in the inferred dust masses. We therefore fix £/ max j = 10 7 in our modeling. 

The limits on adjustable parameters are given in Table [2J The allowed range for U m { n 
is determined by the wavelength coverage of the data used in the fit. For the SINGS galaxy 
sample, it was found that if the photometry extends to A max = 160/im, models with [/ mm > 
0.6 are well-constrained. However, if longer wavelength data are available, we allow the 
possibility of cooler dust, heated by starlight intensities U < 0.6, down to U min = 0.06 if 
Amax = 250/mi, and down to U min = 0.01 if A max > 350/im. 

We observe that for a given set of parameters {(/pah? U m [ n , U mSuX , ol} the model emis- 
sion is multi-linear in {f2*, Md us t, 7}- This allows us to easily calculate the best values of 

M dust , 7} for a given parameter set {(/pah, ^min, f^max, «}• Therefore, when looking for 
the best-fit model in the 7-dimensional model parameter space {f2*, M d , (/pah, 7, ^min ^max, «}, 
we only need to do a search over the 4-dimensional subspace spanned by {(/pah, Knin, t^max, a}- 



12 We try log 10 U max e {3, 4, 5, 6, 7, 8}. 
Table 2: Allowed Ranges for Adjustable Parameters 



Parameter 


min 


max 






Parameter grid used 







n,- 






continuous fit 


M d 





00 






continuous fit 


<?PAH 


0.00 


0.10 






in steps A(/pah = 0.001 


7 


0.0 


1.00 






continuous fit 


^min 


0.7 


30 


when A max 


= 160/im 


unevenly spaced grid a 




0.07 


30 


when A max 


= 250/im 


unevenly spaced grid a 




0.01 


30 


when A max 


= 350/im 


unevenly spaced grid a 




0.01 


30 


when A max 


> 500/im 


unevenly spaced grid a 


O 


1.0 


3.0 






in steps Aa = 0.1 


^max 


10 7 


10 7 






not adjusted 



a The fitting procedure uses pre-calculated spectra for U m \ n G {0.01, 0.015, 
0.01, 0.02, 0.03, 0.05, 0.07, 0.1, 0.15, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 1.0, 1.2, 1.5, 
2.0, 2.5, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 10, 12, 15, 20, 25, 30}. 
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In the case of a fixed £7 max , the search is performed over a 3- dimensional space spanned by 
{qpah, U m i n , a}. In any case, for the computed grid of {qpAR,U TD i a ,U m3X ,a}, the multi- 
dimmensional search for optimal parameters can be performed by brute force, rather than 
needing to rely on a nonlinear minimization algorithm. 

With U max fixed, for each pixel j, the model library is used to search for the model 
parameter vector £ 3 - = M d , ^pah, 7, ^min> ot] that minimizes 

2 _ \- [Fohzj&k)] - ^(model, A fc )] 2 

where -F bs,j(^fc) is the observed flux density and <J\ k j is the 1-cr uncertainty in the measured 
flux density for pixel j at wavelength (see Appendix D for a detailed discussion on how 
tj\ k j is obtained) 

The above procedure yields "best-fit" estimates for the model parameter vector = 
M d , gpAH, 7, ^min, ot) for each pixel j. Each pixel is fitted independently of the remaining 
pixels in the galaxy. Since the final-map pixel size are chosen to Nyquist sample the final- 
map PSF, the camera images are smooth on a pixel scale. The fact that we obtain smooth 
parameter maps is an indication of the stability of the fitting procedure, i.e., even though 
every pixel is modeled independently of its neighbors, the continuity in the images lead 
to continuity in the results. The quoted "best-fit" parameter maps arise from fitting the 
observed flux in each pixels. 



5.1. Properties Derived from the Models 

The infrared luminosity L^j for the dust model in the pixel j is: 

L dJ = P (gpAH,i)M d , i F i , (14) 

where -Po(<7pah) (which depends only weakly on ^pah) is the total power radiated per unit 
dust mass by the model when heated by starlight with intensity U = 1, and Uj is the 



13 In some regions where we do not have complete data coverage or the S/N is low, we may further fix 
some of the parameter values. This situation typically arises in background areas without IRAC coverage, 
in which case we fix f2 + = 

14 We do not use the median maps that could be constructed from modeling adding random noise to the 
observations to estimate the parameter uncertainties 
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mass-weighted mean starlight heating intensity, given by: 



2-ci, 

r/ J - 

^ max 



1 

r/ 3 — u 

'-max w . 



2 — aj 

■U ■ J - 

mm, j 



J 

mm, j 



:i-7i)^«inj+7i< ^ 



ln([/max/f/miii,j) 
ln(C/max/C/mm,j) 



mm,? i_rr . .Ijj 



for a j 7^ 1 , a j 7^ 2 
for CKj = 1 



for aj • = 2 . 



(15) 



Star-forming regions have significant starlight power absorbed by dust grains in regions of 
high starlight intensity, which generally correspond to photodissociation regions (PDRs). 
We will refer to the luminosity radiated by dust in regions with U > Updr as Lpdr, given 
by: 



2- 



-U ■ J 

mm,j 



LpDRj = Po(qPAH,j) M d| j7j X < 



i'n 



'^PDR 



m( f^max / Umin,j 



ln(t/max/t/pDR) 



T 

max 



for a,- 7^ 1 , a,- 7^ 2 



for a,- = 1 



for a. 



(16) 



We take Updr = 10 2 as a plausible cutoff to select dust in high intensity regions (choosing 
another cutoff value would change only the inferred Lpdrj; leaving all the remaining dust 
parameters unaltered). We further define /pdrj as: 



/] 



PDR,j 



^PDRj 



(17) 



The region observed is at a distance D from the observer and ttj is the solid angle of 
pixel j. For each pixel j, the best-fit model vector M d , ^pah, 7, ^min> «}j corresponds to 
a dust mass surface density: 



M, 



Similarly, we can compute the infrared luminosity surface density £i, d j and the surface 
density of dust luminosity from regions with U > Urbr, ^l pdk j, as: 



1 



J L d ,j 



1 



J PDRj 



_D 2 fL 



>PDRj' 



(19) 



The DL07 dust models used here are consistent with the Milky Way ratio of visual 
extinction to H column, Ay/Nn = 5.34 x 10 _22 mag cm 2 /H, for a dust/H ratio SM d /^H^H = 
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0.010 (see Table 3 of iDraine et al.l 120071 ). The dust surface density corresponds to a visual 
extinction (through the disk) 



A v = 0.67 



(lO 5 M kpc~ 



mag 



(20) 



5.2. Global Quantities 

After the resolved (pixel-by-pixel) modeling of the galaxy is performed, we compute a 
set of global quantities by adding or taking weighted means (denoted as (...)) of the quantities 
in each individual pixel of the map. The total dust mass M^ tot , total dust luminosity L^ tot , 
and total dust luminosity radiated by dust in regions with U > Upbr, L^tot, are given by: 

N N N 

M^tot = M d j , L A tot = L d j , L PDR tot = -^pdrj j (21) 

3=1 3=1 3=1 

where the sums extend over all the pixels j that correspond to the target galaxy (see Ap- 
pendix A for the galaxy segmentation procedure). The dust-mass weighted PAH mass frac- 
tion ((/pah), and mean starlight intensity (U), are given by: 

V W , OpAH i Ma 4 y^ y N Uj Ma 

(?pah) ^C^ *; (V) = ^ ' " • (22) 

2^j=l M d,j l^j=l M dtj 

The dust mass-weighted minimum starlight intensity (U m i n ) is given by: 

Uj v _ Ef=l(l-7j)^min,jM d , J - 

Ej-il 1 - lj)Md,j 
The dust-luminosity weighted value of /pdr, (/pdr) is: 

(/pdr) es . (24) 

Alternatively, we can fit the dust model to the global photometry for each galaxyf^l (i.e., a 
single-pixel dust model). The dust parameters obtained from the global-photometry (single- 
pixel) model fit will be compared to the corresponding resolved modeling global quantities 
defined in equations (j2H - EIJ) in £JHJ 



15 We take global photometry to be the photometry within the galaxy mask, i.e., the galaxy regions where 
we can reliably determine the dust luminosity. There is undoubtedly additional emission from material 
outside the galaxy mask, but it cannot be measured reliably. 
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Clearly, the emission summed over the map will not correspond to the emission of dust 
exposed to starlight of the form given by eq. (M- ITOj) since different pixels will have different 
values of f/ m in,j,7j and a.j. For completeness, we could define the dust mass- weighted mass 
fraction heated by a power-law (U > U m \-n ) component (7) as: 

<7> " Z)JLi M iJ ' (25) 

and an "effective power-law exponent" (a), defined as the value that satisfies: 



^PDR.toi 




r r 2-(a> _ TT 2-{a) 
l^max ^PDR 



M") _ ITT . \l-( a ) 



(26) 



Unfortunately, (7) and (a) do not turn out to be useful global quantities, and will, in general, 
differ significantly from the values obtained from the global-photometry (single-pixel) model 
fit. 



5.3. Parameter Uncertainty Estimates 

To estimate uncertainties in the derived dust parameters for each pixel j, we simulate 
data by adding zero-mean random noise 5F(Xk)j, r , for r = 1, 2, 3, N r , to the observed 
flux F(Xk)j in each band, and fit the simulated noisy data. In Appendix [D] we describe 
the statistical construction of the sample F(Xk)j, r of N r random realizations. For the fit 
parameters {a, &,...} G {Q^, M d ,qpAu,l,U min ,a} we have a set of N r + 1 values and we 
calculate the covariance matrix: 

N r 

V ab = V ba = N- 1 J2(a r -a )(b r -b ) , (27) 

r=l 

where do, bo are the best-fit parameter values for the observed fluxes, and the 
best-fit values for the r-th random noise realization. For each model parameter a, the 1-er 
uncertainty is taken to be 

= K 1 / 2 • (28) 

Each random noise realization r = 1,2, ...A^ produces a global quantity via eq. (J2TJ, 
( )22|) . and (|23|) . ( 124)) . and (126]) . We proceed to calculate uncertainties for global quantities us- 
ing equations |27] and |28j In Appendix [E] we describe the details of the parameter uncertainty 
estimation. 
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6. Results 

We construct dust maps with several different angular resolutions. For a given resolu- 
tion, one can only use cameras with the PSF FWHM smaller than the reference PSF. In a 
normal star-forming galaxy, most of the dust mass is at temperatures « 15 — 25 K, with 
vL v peaking at A ~ hc/QkT& « 100 — 160/xm. Reliable estimates of the dust mass therefore 
require long wavelength data, at least out to 160/xm. In the present study we will com- 
pare dust models using maps at the resolution of the PACS160 fFWHM=11.2 "), SPIRE250 



(FWHM=18.2"), and MIPS160 (FWHM=38.8") cameras. As discussed by lAniano et al. 



( 1201 lb the MIPS160 PSF cannot be convolved safely into the SPIRE500 PSF (with FWHM=36.1' 
This implies that MIPS160 is the narrowest PSF that can be used if we want to include 
MIPS160 data into the modeling. A drawback of using the MIPS160 PSF is its extended 
wings, which can cause power radiated by the bright central regions to contribute signifi- 
cantly to the observed surface brightness of the faint outer regions. However, we will see 
in §7 that this effect does not significantly affect dust mass estimates, and discrepancies 
between PACS and MIPS photometry make it important to include the MIPS160 camera in 
the dust modeling. 

In order to generate the H surface density maps (used in the dust mass / H mass ratio 
maps), a value of Aco,2o needs to be chosen. Figure [3] shows the total H maps (first and 
third rows) for both galaxies, for Xqo,2o = 2, 3, and 4, convolved to the MIPS160 PSF. Note 
that the CO maps cover almost all of both galaxy masks, but do not cover the full field of 
view, leading to the box-like step in the H mass maps. 

Using dust maps based on all of the available photometry (with the MIPS160 PSF - 
see Figures St and Eh below), Figure |3] shows the dust/H mass ratios for the different values 
of A C o,2o, for NGC 628 (second row) and NGC 6946 (fourth row). The choice X co ,2o = 4 
gives the smoothest dust/H mass ratios over the galaxies (outside of the central region of 
NGC 6946, which is further discussed in §6.2. ip . As discussed in Sec. I2.3[ we take A COj 2o — 4 
for both NGC 628 and NGC 6946. 
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Fig. 3. — Row 1: Total H surface density for NGC 628 at MIPS160 resolution for X CO ,20 = 2, 3, and 4. 
Row 2: Dust/gas mass ratio derived from the dust map in Figure^, and the gas maps in row 1. Row 3: 
Same as row 1, but for NGC 6946. Row 4: Dust/gas mass ratio derived from the dust map in Figure[9t, and 
the gas maps in row 3. For NGC 628 the dust/gas maps seem best-behaved for Xco,20 = 4 (Figure 3f). For 
NGC 6946 the dust/gas map outside the central region also seem best-behaved for Xco,20 = 4 (Figure 31). 
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Table 3: Galaxy Properties 





NGC 628 


NGC 6946 


Type* 


SAc 


SABcd 


D a (Mpc) 


7.2 


6.8 


Resolution at D 


l"=34.9pc 


l"=33.0pc 


galaxy optical size h (kpc) 


11.0 x 9.9 


11.4 x 9.7 


galaxy mask area A (kpc 2 ) 


274 


307 


galaxy mask radius \J A/ir (kpc) 


9.4 


9.9 


M(H I) c (M Q ), total in map 


(3.7 ±0.2) x 10 9 


(5.5 ±0.3) x 10 9 


M(H 2 ) d (Af©), total in map 


(2.5 ±0.6) x 10 9 


(9.7 ±2.7) x 10 9 


M(H I) (M 3 ), within mask 


(1.5 ±0.1) x 10 9 


(2.1 ±0.1) x 10 9 


M(H 2 ) d (M ), within mask 


(2.0 ±0.2) x 10 9 


(8.6 ±0.5) x 10 9 


Mh (Mq), within mask 


(3.5 ±0.3) x 10° 


(10.7 ±0.6) x 10 9 


log 10 (O/H) + 12 e 


8.27±0.12 


8.37 ±0.06 


L(Ha) f (L Q ), total 


1.88 x 10 7 


1.00 x 10 8 


SFRf (Mgyr- 1 ), total 


0.7 ±0.2 


4.5 ± 1.2 


F„(IRAC3.6) (Jy), within mask 


0.83 ±0.10 


3.13 ±0.52 


F y (IRAC4.5) (Jy), within mask 


0.57 ±0.06 


2.22 ±0.32 


i^(IRAC5.8) (Jy), within mask 


0.98 ±0.23 


4.9 ± 1.2 


F^IRACS.O) (Jy), within mask 


2.61 ±0.44 


13.2 ±2.4 


F„(MIPS24) (Jy), within mask 


3.04 ±0.33 


18.9 ±2.0 


i^(MIPS70) (Jy), within mask 


31.5 ±9.8 


196 ±55 


F^MIPSIGO) (Jy), within mask 


104 ± 20 


420 ± 127 


F„(PACSS70) (Jy), within mask 


39 ± 12 


245 ± 58 


F„(PACSS100) (Jy), within mask 


73 ±21 


433 ± 122 


F^PACSSIGO) (Jy), within mask 


110 ±20 


526 ± 126 


F I/ (SPIRE250) (Jy), within mask 


59.5 ±7.0 


247 ± 29 


F„(SPIRE350) (Jy), within mask 


27.4 ±3.3 


101 ± 12 


F„(SPIRE500) (Jy), within mask 


10.5 ± 1.4 


35.4 ±4.2 


Md (Mq), within mask 


(2.9 ±0.4) x 10 7 


(6.7 ±0.6) x 10 7 


100 x M d /M H , within mask 


0.82 ±0.17 


0.63 ±0.09 


Ld (Lq), within mask 


(6.8 ±0.4) x 10 9 


(3.3 ±0.3) x 10 10 


Lpdr (Lq), within mask 


(7.9 ± 1.4) x 10 8 


(4.7 ±0.8) x 10 9 


(U m i n ), within mask 


1.6 ±0.3 


3.1 ±0.6 


U, within mask 


1.7±0.3 


3.6 ±0.6 


(/pdr), within mask 


(11.6 ± 1.3)% 


(14.2 ±2.8)% 


(<7pah), within mask 


(3.7 ±0.3)% 


(3.6 ±0.4)% 



a From Kennicutt et al. ( 20111 . 

b Major and minor radii at ^s=25 mag arcsec -2 isophote, from Kennicutt et al. (2011) 
c We use the "Natural" (NA) weighting maps fsee Walter et al. (2008) for details). 
d From Lerov et al. (2009), for A co ,20 = 4 (see text). 

e PT05 HII region abundances, unweighted average, from Moustakas et al. (2010) . 
f From Kennicutt et al. f 20081 . uncorrected for internal extinction. 
9 From Calzetti et al. ( 20101 . based on Ha+24/im. 
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6.1. NGC 628 

6.1.1. Maps of Gas and Dust 

NGC 628 (= M74), at a distance D = 7.2 Mpc, is classified as SAc. With major and 
minor optical diameters of 10.5 and 9.5 arcmin, it is well resolved even by the MIPS160 
camera. 



The star formation rate is estimated to be 0.7 ± 0.2M Q yr 1 f Calzetti et alJlioiol ). Two 



supernovae have been observed in NGC 628: SN 2002ap (Type Ic) and SN 2003gd (Type 
II-P). 



The total H I mass from 21-cm observations is M(H I) = (3.7±0.2) x 1O 9 M flWalter et al. 



20081 ) and the tota 



H 2 mass estimate d from observations of CO 2-1 is M(H 2 ) = (2.5 ±0.6) x 



lO 9 (Xco,2o/4)M (jLeroy et al.ll2009l ). The fact that the adopted value of Xco is larger than 
the value Xco, 20 ~ 2 found for resolved CO clouds in the Milky Way will be discussed in 
Section [10] below. 

Figure H] shows the resulting maps of H surface density (panels a-c), dust surface den- 
sity £j\/ d (panels d-f), dust luminosity surface density Si d (panels g-i), and dust/H mass 
ratio (panels j-1), obtained by fitting photometry with PACS160, SPIRE250, and MIPS160 
resolution, using all the compatible cameras in each case. The dust models with PACS160 
resolution show clear spiral structure, but the noise in the smaller pixels is such that the 
dust is only reliably detected at surface densities S A/d > 10 5,0 M Q kpc -2 , corresponding to 
Av ^ 0.7 mag. If the mapping is done at the resolution of the SPIRE250 camera, the dust 
is reliably detected for Em,, > 1O 4,5 M kpc -2 , corresponding to Ay > 0.2 mag; at this reso- 
lution, the spiral structure is still visible. Maps made at MIPS160 resolution are the most 
sensitive, because of the larger pixel size, and the fact that they use data from all of the 
cameras, including the MIPS160 camera. Unfortunately, the lower resolution of these maps 
(36" FWHM) largely washes out the spiral structure which is visible in higher resolution 
maps of NGC 628. Nevertheless, imaging at MIPS160 resolution allows reliable detection of 
dust at surface densities as low as 1O 4 ' 3 M kpc -2 , or A v > 0.14 mag. 

We note that the dust/H mass ratio maps change significantly when the modeling is done 
at the different resolution/ camera combination (see last row of Fig. H]). These discrepancies 
are mainly due to the low sensitivity of the PACS cameras in the low-surface brightness 
areas (compare the MIPS and PACS flux in the outer parts of the galaxies in Figures [1] and 
[T]), and discrepancies in the high surface brightness areas. In §7 we discuss the PACS-MIPS 
photometry discrepancies further. 
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Fig. 4. — NGC 628 at the resolution of PACS160 (left), SPIRE250 (center), and MIPS160 (right). PACS160 
resolution models are based only on IRAC, MIPS24, and PACS data; SPIRE250 resolution models are based 
on IRAC, MIPS24, PACS and SPIRE250 data; MIPS160 resolution models are based on all (IRAC, MIPS, 
PACS, SPIRE) data. Row 1: surface density of H (both HI and H 2 ) for A co ,20 = 4 (sec text). Row 2: 
estimated dust surface density Sj\/ d (see text). Row 3: dust luminosity surface density £z, d . Row 4: dust/H 
mass ratio over the main galaxy. The irregular white contour is the boundary of the "galaxy mask" (see 
text). White circles are selected apertures (see Fig. [5]). 
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Fig. 5. — NGC 628 dust models at the resolution of PACS160 (left), SPIRE250 (center), and MIPS160 
(right). Top row: PAH abundance parameter gpAH- Middle row: diffuse starlight intensity parameter C/ m i n . 
Bottom row: PDR fraction /pdr. 



6.1.2. Maps o/gpAH an d Starlight Parameters 

Figure [5] shows maps of the PAH abundance <2pah (panels a-c), the starlight intensity 
parameter U m i n (panels d-f), and the PDR fraction /pdr (panels g-i), over the "galaxy mask" 
region where the galaxy is well-detected. The galaxy mask for NGC 628 has a diameter of 
~0.16°, or 20 kpc © 7.2 Mpc. 
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The PAH abundance parameter qpah, shown in Figure [5^-c is remarkably uniform over 
the region where it can be reliably estimated. In Figure |5b, (/pah varies from a high of ~ 0.05 
a few kpc from the center down to ~ 0.035 near the edge of the galaxy mask. If there is a 
radial gradient in (/pah? it is we ak, consistent with the weak gradients found for the SINGS 



sample (including NGC 628) bv iMunoz-Mateos et all (120091 ). 



At PACS160 resolution, the starlight intensity parameter U min (see Fig. [5H) varies be- 
tween ~0.6 and ~3 over most of NGC 628, following the galaxy structures, but near the 
edges of the galaxy mask U min appears to rise. This is because the reduced signal/noise 
results in PACS70/PACS160 ratios that appear to be anomalously high, leading to high 
inferred values of U m i n for some pixels. This is probably the result of low PACS160 fluxes for 
those pixels, making it appear that the dust is rather warm. We see that when SPIRE250 
data is introduced (Fig. [5fe), we have many fewer high values of t/ m ; n near the edge of the 
galaxy mask, and the ?7 m i n values in the brighter regions appear well-behaved. This continues 
when the MIPS160, SPIRE350, and SPIRE500 data are brought into the fit in the MIPS160 
resolution image (Fig. [5f). 

The bottom row of Figure E] (panels g-i) shows /pdr, the fraction of the dust luminosity 
coming from dust heated by starlight intensities U > 100, which we expect to be associated 
with star-forming regions. At PACS160 resolution, most of the galaxy has /pdr ~ 0.03 and 
some small bright regions (for example the center of aperture 2 in Figure |8]) have higher 
values, up to /pdr ~ 0.3. These regions with high values of /pdr generally coincide with 
Ha peaks, and also with the regions of the highest dust luminosity per area (see Fig. Hg). 
As we shift to coarser modeling (i.e, using PSFs with larger FWHM), the /pdr peaks are 
smoothed, and the general galaxy pixels tend to have larger /pdr values: the dynamic range 
of values of /pdr decreases. At MIPS160 resolution, most of the galaxy has /pdr > 0.05. 
The overall (dust luminosity-weighted) mean for the galaxy is (/pdr) = 0.116: 11.6% of the 
total IR power is radiated by dust in regions where U > 100. 

Unfortunately, at MIPS160 resolution the arm structure of both galaxies is not clearly 
resolved (see the MIPS 160 images of the galaxies in panel (f) of Figures [1] and [2]) . Therefore 
we cannot reliably study variation of the dust parameters between arm and interarm regions. 
A subsequent work (Hunt et al. 2012, in prep.) will study radial variations in the model 
parameters. 
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Fig. 6. — Ratio of model/observed intensity at A=250/im (top row), A=350/im (middle row), and 
A=500 / um (bottom row) for NGC 628. Left column: PACS160 PSF, model constrained only by IRAC, 
MIPS24 and PACS. Center column: SPIRE250 PSF, model constrained by IRAC, MIPS24, PACS, and 
SPIRE250. Right column: MIPS160 PSF, model constrained by all 13 cameras. The model in the center 
column does a fairly good job in predicting the emission at A=350^m , and A=500/im , except near the edge 
where the S/N is low. The model in the right column is in excellent agreement with all three SPIRE bands. 
In all panels the model predicted and observed intensities have been convolved to a common PSF before 
taking the ratios. 
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6.1.3. Comparison Between Observed and Modeled Flux Densities 

How well does the dust model reproduce the SPIRE photometry? Figure [6] shows 
the ratios of model-predicted intensity to observed intensity at A=250, 350 and 500/imfor 
dust models obtained by fitting photometry with PACS160, SPIRE250, and MIPS160 res- 
olution. In order to make the comparison, we degrade either the observed image or the 
model-predicted image to a common resolution (i.e., when the modeling is done at PACS160 
resolution, we convolve the model-predicted SPIRE250 image to the SPIRE250 PSF, and 
when the modeling is done at MIPS160 resolution, we convolve the observed SPIRE250 im- 
age to the MIPS160 PSF). Except near the edge of the galaxy mask (where the low S/N in 
the PACS data becomes an issue), the modeling tends to owerpredict the SPIRE500 pho- 
tometry - there is no evidence for a significant mass of very cold dust radiating only at the 
longer wavelengths. 

When we attempt to model at PACS160 resolution (using only IRAC, MIPS24, and 
PACS data to constrain the model), the model predictions at SPIRE250, SPIRE350, and 
SPIRE500 do not agree very well with observations (see Fig. [6^, d, and g). If we coarsen 
the modeling to SPIRE250 resolution and add SPIRE250 to the model constraints, we now 
reproduce the SPIRE250 image (not surprising) and the model predictions at 350/imand 
500/imare within ±30% and ±50% of the SPIRE observations respectively over most of the 
galaxy mask. If we further coarsen the modeling to the MIPS160 PSF, and use all the data 
to constrain the model, we find good agreement at all SPIRE bands (see Fig. [6b, f, and i), 
with only SPIRE500 slightly over-predicted by ~ 10% in some regions. 

Figure Et, f, and i show that when all of the cameras are used to constrain the modeling 
(and with the improved S/N of the larger MIPS 160 pixels), the model emission is gener- 
ally in good agreement with the SPIRE imaging - for all three SPIRE bands the ratio of 
model/observation is close to 1 over most of the galaxy, with significant departures only just 
at the edge of the galaxy mask, where observational errors are likely to be the explanation. 

6.I.4. Global SED Fitting 

Figure [7| shows global SEDs for NGC 628. The observed photometry is represented by 
rectangular boxes (Spitzer (IRAC, MIPS) in red; Herschel (PACS(S), SPIRE) in blue). The 
model convolved with camera response function (i.e., expected camera photometry for the 
model) is represented by diamonds (Spitzer in red, Herschel in blue). The two rows show 
the global SED for NGC 628. The observed global photometry for NGC 628 is the same 
across the two top rows, but the models differ. Each black curve is a fitted model spectrum. 
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Fig. 7. — Model SEDs for NGC 628. Black line: total model spectra. Cyan line: stellar contribution. Dark 
red line: emission from dust heated by the power-law U distribution. Dark green line: emission from dust 
dust heated by U = U min - Solid-line rectangles: observations used in the fit (red: Spitzer (IRAC, MIPS): 
blue: Herschel (PACS(S), SPIRE)). Dashed-line rectangles: observations not used in the fit (same color 
scheme as solid- line rectangles). Diamonds: model convolved with camera response function (i.e., expected 
camera photometry of the model) (red: Spitzer; blue: Herschel). Top row: Global SED compared with single- 
pixel models. Second row: Global SED, for multi-pixel model. In the left column, only IRAC, MIPS24, 
and PACS are used in the fit (MIPS70,160 and SPIRE are not used), yet the model nevertheless falls close 
to the observed SPIRE fluxes. In the center column, IRAC, MIPS24, PACS, and SPIRE250 are used (i.e., 
MIPS70,160 and SPIRE350, 500 are not used), and the agreement with SPIRE is improved. In the right 
column, all IRAC, MIPS, PACS, and SPIRE data are used to constrain the model, and the agreement with 
all SPIRE bands is excellent, although the model slightly overpredicts the emission at 350/miand 500/xm. 
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The diamonds show the fitted model spectrum convolved with the instrumental response 
function for each camera, i.e, the expected photometry for the fitted model. 

The top row shows "single pixel" models for the galaxy based on (a) IRAC, MIPS24, 
and PACS data only, (b) IRAC, MIPS24, PACS, and SPIRE250, and (c) IRAC, MIPS, PACS 
and SPIRE (all the cameras). 

The second row shows the predicted model SED obtained by fitting a dust model to 
each resolved pixel, and then summing the emission over all the pixels. These predicted 
SEDs differ from the previous "single pixel" predicted SED because the dust modeling is 
a non-linear process. The total dust mass predicted in this way is, however, similar to the 
"single pixel" prediction (see §8 for details). When the modeling is done at lower resolutions 
(MIPS160) the "single pixel" and map-averaged quantities are in closer agreement. 

6.1.5. Fitting in Selected Apertures 

Figure [8] shows SEDs for four circular apertures located on NGC 628, sampling a wide 
range of surface brightnesses, S id « (0.3—9) x 1O 7 L kpc -2 (see Fig.HJfor aperture locations). 
Aperture 4 is located partially outside the galaxy mask, where the single pixel dust modeling 
is not reliable, but the improved S/N of a large aperture allows reliable determination of the 
dust model parameters. We fit the DL07 model to the summed flux within each aperture. 
The top row shows the SED for a 60" (diameter) circular aperture centered on the galaxy 
nucleus. The second row shows the SED for a 60" aperture located on a bright spot on 
the spiral arms. We observe that the PACS photometry is larger than the corresponding 
MIPS photometry in the high surface brightness apertures 1 and 2. When both data sets 
are included in the modeling, the dust model gives values intermediate between PACS and 
MIPS. The third and fourth rows show SEDs in 80" and 120" apertures further from the 
center. We employ larger apertures in order to obtain reasonable S/N in these fainter regions. 
In aperture 4, the MIPS and PACS photometry differ significantly. When SPIRE, MIPS70, 
and MIPS160 are not used (e.g., in the PACS160 resolution modeling, Fig. [Hj) the high 
PACS70/PACS160 ratio causes the model to infer very high values of £/ m i n and low values 
of EM d , and hence to underpredict the SPIRE photometry. When we include SPIRE250 
(Fig. [Hk) the model can reproduce SPIRE250, but continues to underpredict SPIRE350 
and SPIRE500. Finally, when SPIRE350,500 and MIPS70,160 are added as constraints 
(Fig. [H), the modeling improves dramatically in aperture 4 and other low-brightness areas, 
reproducing most of the data, and making it clear that PACS70 is an outlier. 

The astute reader will note that the estimated uncertainties for, e.g., the PACS 100 
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Fig. 8. — Model SEDs for four selected apertures on NGC 628. Sec Fig. [4] for aperture location, and 
Fig. [7] for an explanation of the color coding. Top row: 60" (circular) aperture centered on the galaxy 
nucleus. Second row: 60" aperture located on a bright spot on the spiral arms. Third row: 80" aperture in 
a mid- luminosity region. Bottom row: 120" aperture in a low-luminosity region. 
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global photometry differ between the columns (i.e., for fluxes extracted after convolving to 
different PSFs). For each PSF, we estimate the noise per pixel based on the pixel statistics 
in the background region (see Appendix [D]). We then make a simple assumption concerning 
the pixel-to-pixel noise correlation. The fact that the uncertainty estimates for the aperture 
fluxes depend on the PSF is an indication that our assumption about the correlated compo- 
nent of the noise is imperfect. This is only an issue for the faint, low S/N data, such as the 
PACS fluxes in Apertures 3 and 4. 

In aperture 4 (Fig. [S]last row), the model uses 7 = 1, i.e., the dust is heated almost 
entirely by a power law U distribution, but with very high values of 7. This corresponds to 
a broad distribution of starlight intensities in the U > t/ m i n range, with very little power in 
the high intensity range (PDR). This may be an artifact arising from the large photometric 
uncertainties. 



6.2. NGC 6946 



NGC6946, an active star-forming galaxy, is classified as SA Bcd. At a d i stanc e D = 
6.8 Mpc, the total H I mass is Af(HI) = (5.5 ± 0.3) x 1O 9 M JWalter et al.l l2008h . with 
42% of the H I falling within the galaxy mask. The IE m ass within the galaxy mask is 
M(H 2 ) = (8.6 ±0.6) x lO 9 (X co ,2o/4)M JLerov et allkopgh . As discussed previously, we 
adopt Xcopn = 4 as a global estimate. The star formation rate is estimated to be 4.5M yr _1 



( jCalzetti et al.l 120101). NGC 694 6 is remarkable for hosting at least 9 supernovae over the 
past century ( iPrieto et al.l 120081 1 . 



The variable carbon star V0778 Cyg, located near RA=309.044, Dec= 60.082 (slightly 
off the bottom left corner of the maps shown in figures 191 - [TT1) . saturates the IRAC and MIPS 
detectors. The associated image artifacts affect the background estimation and subtraction, 
making the modeling less reliable in the bottom left corner of the maps. 



6.2.1. Dust and H Mass Maps 



Figure M shows maps of the gas in NGC 6946 obt ained from the TH INGS 21-cm map 
f lWalter et aDl2008h and the HERACLES CO 2-1 map JLerov et al.ll2009h . convolved to the 
resolution of PACS160, SPIRE250, and MIPS160. The second row shows dust maps obtained 
by the present study. The dust map obtained using the PACS 160 PSF has dust clearly visible 
only out to a radius of ~250" (8kpc @ 6.8 Mpc), but the 11" (500 pc) FWHM of the beam 
resolves the spiral structure. The arm/interarm contrast in dust surface density appears to 
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be approximately a factor of ~2 in this image. Introducing SPIRE 250/im data requires the 
PSF to be broadened, making the spiral arms less apparent, but allows dust to be detected 
out to larger radii, notably in the northern spiral arm. The MIPS160 resolution image clearly 
shows dust in the northern spiral arm out to a distance of ~15kpc from the nucleus. 

The dust luminosity/area Sx, d (Figure [U bottom row) peaks at ~10 101 L Q kpc~ 2 in the 
nucleus (see the PACS160 resolution map), and can be followed down to 1O 7 O L kpc~ 2 

in the MIPS160 resolution map. Similarly, the IR luminosity/area contributed by PDRs 
peaks at £l pdr ~ 10 9 L Q kpc -2 near the center, and is reliably measured down to ~ 1O 5,6 L kpc~ 

The dust/H mass ratio ratio shown in Figure[9j-1, calculated assuming Aco,2o = 4, shows 
a pronounced minimum of ~ 0.005 in the central kpc. In this region the gas is primarily 
molecular, and the estimated dust/H mass ratio is therefore sensitive to the value assumed 
for X co . NGC6946 has approximately solar metallicity, and we expect the dust/H mass 
ratio to be ~0.01. The surprisingly low dust/H mass ratios found in the center of NGC6946 
indicate that the value of Xco near the center should be about a factor of ~ 2 smaller than 
the value X CO) 2o = 4 adopted for the bulk of the galaxy. This conclusion is consistent with 
i nterferometric studies of GM Cs which indicate Xqo,2o ~ 1-25 near the center of NGC6946 
f lDonovan Meyer et al.ll2012l ). using virial mass estimates for individual GMCs. 



6.2.2. Maps of qpAH and Starlight Parameters 

The PAH abundance ^pah in NGC 6946 rises from a minimum of < 1% near the nucleus, 
reaching levels of ~(4 ± 1)% at distances ~l-8kpc from the center, ultimately appearing 
to decline at the edge of the detectable region. The central minimum in (?pah is real (it 
is independent of uncertainties in the appropriate value of Xqo). The local minimum of 
gpAH at the center may be the result of destruction of PAHs by processes associated with 
star formation (there is no evidence of AGN activity in NGC 6946). We note that there 
are a number of other local minima of (/pah evident in Figure [TUh-c; just as for the central 
minimum, these extranuclear minima tend to coincide with peaks in dust luminosity surface 
brightness (see Figure [9g,h) and peaks in /pdr (see Figure [10],k,l). These are both 
indicators of luminous star-forming regions, which can be expected to coincid e with H II 



regions around hot stars, which appear to destroy PAHs (e.g.. |Povich et al.ll2007l ). Supernova 
blastwaves are presumed to also destroy PAHs. 

Because the outer falloff in gpAH is occurring as the signal/noise ratio is falling to low 
values, it is not certain whether the observed decline is real or an artifact of different levels of 
background subtraction at different wavelengths. However, the decline in (/pah in the outer 
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regions persists in the MIPS160 resolution map. These variations in qpah will be examined 
in more detail in future studies. 

The diffuse starlight intensity parameter U m - m (Fig. [TO], middle row) shows a general 
decline with increasing distance from the center. The PACS160 resolution map of U m \ n is 
somewhat noisy, particularly in the outer regions, but the MIPS 160 resolution map of U mm 
shows a systematic decline from values as high as U m \n ~ 10 in the central ~ 500 pc down 
to U min ~ 0.15 in the outer regions, ~ lOkpc from the center. Even lower values of U min are 
estimated for some pixels, but this is seen only at the lowest signal/noise levels. 

The bottom row of Figure [TOl shows /pdr- The behavior of /pdr in NGC 6946 is similar 
to NGC 628. Bright complexes coincide with local maxima of /pdr- Outside very bright 
complexes /pdr is quite smooth across the galaxy. 

6.2.3. Comparison Between Observed and Modeled Flux Densities 

The model-predicted/observed flux ratios for NGC 6946 (shown in Figure [TT]) behave 
similarly to those of NGC 628 (see Figure [6]). When long wavelength data is not used in the 
dust model fit, the modeling tends to overpredict the intensity at the longer wavelengths. 
However, when the SPIRE photometry is included in the model constraints, the overpre- 
diction is greatly reduced. In the case of MIPS160 PSF (all bands used), the SPIRE250 
photometry is reproduced by the model within 10%, and SPIRE500 tends to be overpre- 
dicted by less than 15%. 

6.2.4. Global SED Fitting 

Figure [T2] shows the global SED for the galaxy NGC 6946 (similar to Figure [7] for 
NGC 628). The top row shows the SED for a "single pixel" model and the second row shows 
the predicted model SED obtained by summing over the model SED for individual pixels. 
Again, we observe small differences between the global modeling (top row) and resolved 
studies (bottom row) due to the non-linear behavior of the modeling process. It is seen 
that when the modeling is done using the MIPS 160 PSF, the dust mass estimated from the 
multi-pixel model is in good agreement with that obtained for a single-pixel model. 



-41 - 



6. 2. 5. Fitting in Selected Apertures 

Figure [13] shows SEDs for four circular apertures located on NGC 6946, sampling a 
wide range of surface brightnesses, « (0.07— 13) x 10 8 L Q kpc -2 , similar to Figure [8] for 
NGC 628. The top row shows the SED for a 60" diameter circular aperture centered on the 
galaxy nucleus. The second row shows the SED for a 60" aperture located on a bright spot 
on the spiral arms. The third and fourth rows show SEDs in 120" apertures further from the 
center. Aperture 4 is located completely outside the galaxy mask, but the S/N of a large 
aperture allows a reliable determination of the dust model parameters. 

The nuclear emission is strongly peaked. The Herschel photometry for the central 
aperture noticeably decreases when the image is convolved to MIPS 160 resolution, and the 
estimated dust luminosity surface brightness in the 60" extraction aperture drops from 2.2 x 
10 9 L Q kpc- 2 (Figs. H3t,b) to 1.3 x lO 9 L kp C - 2 (Fig. US). 

We observe that the PACS photometry in the high surface-brightness nucleus tends to 
generally be larger than the corresponding MIPS photometry - compare the PACS and MIPS 
fluxes in Figure [T3b. In apertures 2, 3, and 4, the PACS and MIPS photometry (at common 
resolution - see Figures [TBT . i, and 1), the PACS and MIPS photometry is in fairly good 
agreement. The disagreement in the PACS vs. MIPS photometry is larger for NGC 6946 
than for NGC 628 (see Table [6] below) , presumably related to the much higher peak surface 
brightnesses present in NGC 6946. At MIPS160 resolution, = 9 x lO 7 L kpc -2 in 
NGC 628 aperture 1 is only 7% of S Ld in NGC 6946 aperture 1. 



-42 - 



PACS160 PSF 



SPIRE250 PSF 



MIPS160 PSF 



60.3 4iGC'f>9~:& 




60.3 4JGG6!"Kfi- 



60.3 4JGC6946 100 x Dust Mass / H Mass 




llllIlL 



■3.36 " 

lapa - 
lg.se _ 
I2.00 

■l.68~ 

-1.41 " 
■ 1.19 - 

1'i.bo 
0.04 

'o.'?i 

■0.59- 
■0.50 _ 
lo.42 
lo.35 

loio " 

lo.25 - 



60 - 1 
i J 1 

309 




rOO- x- Du-5t- Mas:s /■■ H -Mass- 




— 1 



308.8 308.6 
RA (deg) 



308.8 308 6 

RA (deg) 



i09 308.8 308.6 

RA (dog) 



Fig. 9.— NGC 6946 at the resolution of PACS160 (left), SPIRE250 (center), and MIPS160 (right). Row 
1: surface density of HI and H2 for Xco,20 = 3 (see text). Row 2: dust surface density Smj- Row 3: dust 
luminosity surface density Si d . Row 4: dust / H ratio over the main galaxy. The irregular white contour 
is the boundary of the "galaxy mask", within which the signal/noise ratio is high enough to obtain reliable 
estimates of dust and starlight parameters. White circles are selected apertures (see Fig. [T3|) . 
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Fig. 10. — NGC 6946 dust models at the resolution of PACS160 (left), SPIRE250 (center column), and 
MIPS160 (right). Top row: PAH abundance parameter <7pah- Middle row: diffuse starlight intensity param- 
eter [/ m i n . Bottom row: PDR fraction /pdr- 
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Fig. 11. — Similar to Fig.H but for NGC 6946. Ratio of model intensity / observed intensity for NGC 6946 
at A=250/im(top row), A=350/im (middle row), and A=500^m (bottom row). Left column: PACS160 PSF, 
model constrained only by IRAC, MIPS24 and PACS. Center column: SPIRE250 PSF, model constrained 
by IRAC, MIPS24, PACS, and SPIRE250. Right column: MIPS160 PSF, model constrained by all 13 
cameras. The model in the center column does a fairly good job in predicting the emission at A=350/im , 
and A=500/mi, except near the edge where the S/N is low. The model in the right column is in excellent 
agreement with all three SPIRE bands. 
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Fig. 13. — Model SEDs for four selected apertures on NGC 6946. See Fig. [9] for aperture location, and 
Fig. [7] for an explanation of the color coding. Top row: 60" (circular) aperture centered on the galaxy 
nucleus. Second row: 60" aperture located on a bright spot in a spiral arms. Third row: 120" aperture in a 
mid-luminosity region. Bottom row: 120" aperture in a low-luminosity region. 
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7. Importance of SPIRE Photometry 

We investigate the effect of only using certain sets of cameras in our resolved (multi- 
pixel) modeling, since this situation will arise in galaxies that were only observed with some 
of the instruments. Dust estimates for NGC 628 and NGC 6946 based on different camera 
combinations are shown in Figure HU where we compare the dust mass estimate with our 
"gold standard", the dust mass estimate using all cameras, at MIPS160 resolution^!. 



M Du „ estimate for different cameras and PSFs 




PSF and cameras used 

Fig. 14. — Total dust mass (summed over all pixels) estimated for NGC 628 and NGC 6946 using 
different camera/PSF combinations (see text). Masses are given relative to the "best estimate" Md obtained 
summed over all pixels using all wavelengths and MIPS160 resolution. 



In Figures 1X411171 the horizontal axis represents different combinations of PSF and cam- 
eras. The nomenclature convention is as follows: The first characters describe the PSF 
used; P160, S250, S350, S500, M70, and M160 refer to PACS160, SPIRE250, SPIRE350, 
SPIRE500, MIPS70, and MIPS160 respectively. 

The first group of 3 digits refers to the MIPS24, MIPS70, and MIPS160 cameras, with 
1 indicating usage. The next set of three digits similarly indicates usage of the three PACS 
bands, and the final three digits correspond to usage of the SPIRE cameras. For example, 
S350 110 111 110 uses MIPS24, MIPS70, PACS 70, PACS100, PACS 160, SPIRE250, and 
SPIRE350 cameras at SPIRE350 PSF. We always use the Scanamorphos data reduction for 
PACS. 



16 Even though M1PS160 PSF is our lowest resolution modeling, PACS - MIPS discrepancies made it 
important to be able to retain the MIPS160 camera in the modeling, and, as shown in §8, using lower 
angular resolution does not introduce a major bias in the modeling. 



-48- 



The ten right-most columns correspond to the MIPS160 PSF, and the remaining eight 
left- most columns have different PSFs (sorted by ascending FWHM left to right). The mod- 
eling done at MIPS160 resolution uses different camera combinations, allowing one to exam- 
ine, at fixed angular resolution, the effects of different wavelength coverage and PACS/MIPS 
discrepancies. 

Figure [14] shows that the dust mass estimates for these two galaxies appear to be quite 
good when we use the three SPIRE bands: the total dust mass estimated without using 
MIPS160 (at MIPS160 or SPIRE500 PSF) is off by only ~10%. It therefore appears that re- 
solved maps of dust in nearby galaxies can be reliably obtained at the resolution of SPIRE500. 
The dust mass estimates appear to be good provided at least two SPIRE bands (SPIRE250 
and SPIRE350) are employed: the total dust mass estimated without using SPIRE500 and 
MIPS160 is off by 30% for NGC628, and 22% for NGC6946. The S250 100 111 100 dust 
map yields a global dust mass that exceeds the best estimate by ~38%. For some pur- 
poses this ~38% loss of accuracy will be acceptable, as it makes possible a dust map with 
FWHM=18.2". 

In order to map the dust at higher angular resolution (PACS160 PSF), we can only use a 
subset of all the cameras available (IRAC, MIPS24, PACS). Dust maps made at the resolution 
of PACS160 are less reliable, for three reasons: (1) with smaller pixels, the signal/noise of the 
image is poor in the faint regions, (2) there seem to be significant systematic uncertainties 
associated with the PACS photometry, (3) the dust models are unconstrained at A > 160/im, 
and (4) our PACS maps are less sensitive than the MIPS counterparts. The non-linearity of 
the dust model fitting procedure can result in large errors in dust mass in the low S/N regions. 
Dust masses estimated using maps at PACS160 resolution (using IRAC, MIPS24, and PACS 
data only) can overestimate the dust masses by factors ~ 1.4 — 1.8. If greater accuracy is 
required, PACS 160 resolution should only be used in the highest surface brightness regions, 
where signal/noise issues should be less critical. 

It has sometimes been argued that dust masses estimated using only data at wave- 
lengths A < 160/im (e.g., those based on MIPS photometry only) are suspect because the 
observations are insensitive to "cold" dust at t emperatures Th S, 1 2 K. The SINGS sample 



included 17 galaxies with SCUBA 850/mi, and lDraine et al.l (120071 ) found that reliable dust 
mass estimates could be obtained with the DL07 dust model without using data longward 
of 160/xm : these galaxies did not appear to contain significant masses of "cold dust" . The 
SPIRE photometry for NGC 628 and NGC 6946 confirms this: these galaxies do not appear 
to contain substantial masses of cold dust. In fact, we find that models fitted to the PACS 
photometry alone tend to predict more emission at long wavelengths than is observed. This 
is evident in Figures [6^ and[6]i, and Figures [TTa andCETH, which show that DL07 dust models 
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based only on A < 160/im data (IRAC, MIPS24, and PACS) tend to overpredict the actual 
fluxes at 250 and 500/imin the central 6 arcmin (12kpc). This is also seen in Figures [7U and 
IT2H . which compare the observed global SED with the sum of the SEDs from a multipixel 
model for each galaxy. The model (the diamonds) overpredicts the SPIRE data by factors 
of ~1.1, 1.2, and 1.2 at 250, 350, and 500/im. This is rather good agreement, and indicates 
that there is no substantial amount of very cold dust present in either galaxy. In fact, when 
the SPIRE data are used to constrain the fit, the best-fit models actually raise the starlight 
intensities slightly in order to match both the PACS, MIPS, and SPIRE photometry, reduc- 
ing the dust mass estimate. Note the very good agreement between observed and model 
SEDs in Figures [7f and [T2T . where all data are used to constrain the multipixel fit. 

The main difference in the estimated dust masses in the different PSF/cameras com- 
bination is due to the different values of U m \ n found when we do not use all the cameras 
available. Figure [151 shows (U m i n ) (defined in eq. [23]) , the dust-mass weighted value of the 
starlight intensity parameter U m i n relative to the best estimate for (U m i n ) - the estimate 
obtained using all the data at MIPS 160 resolution. 



< U m in> estimate for different cameras and PSFs 




Fig. 15. — (t/min) (weighted mean over the pixels) estimated using different cameras and PSFs relative to 
the best estimate (obtained with all cameras and MIPS160 PSF). 



8. Maps v.s. Global Photometry 



The KINGFISH sample consists of galaxies that are near enough to be well-resolved, 
enabling us to make maps of dust and starlight properties. However, because the dust 
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modeling used here is a nonlinear procedure, the total dust mass estimate, for example, will 
depend on what resolution is employed. In the preceding sections, we usually made dust 
maps using the best angular resolution possible with each chosen camera set used, with the 
limiting resolution therefore determined by the longest wavelength data used in the fit (or 
the inclusion of MIPS160). 

Herschel is frequently used to observe distant galaxies that are unresolved at the longest 
wavelengths. For the KINGFISH galaxies, we will therefore also estimate dust masses using 
the global photometry (the flux from within the galaxy mask), to see how dust mass estimates 
based on global photometry compare with the best estimates for a well-resolved galaxy. 



1.7 




PSF and cameras used 

Fig. 16. — Dust mass estimates based on global photometry for NGC 628 and NGC 6946, using different 
camera combinations, relative to the sum of the dust mass (obtained from pixel-by-pixel modeling of the 
resolved galaxy) in the MIPS160 PSF using all the cameras (i.e., relative to the best total dust estimate for 
the galaxies). In principle estimates based on global photometry should not depend on the chosen PSF, only 
on the camera combination. In practice, flux scattered out of the galaxy mask (into the PSF's extended 
wings) changes the global photometry estimates, leading to small differences in the parameter estimates 
(e.g., sec differences between P160 100 111 000 and M160 100 111 000 modeling). 



Figure [16] shows the dust mass estimated using the DL07 dust models to fit the global 
SED, divided by the "gold standard", i.e., the sum of the dust mass obtained from pixel-by- 
pixel modeling of the resolved galaxy in the MIPS 160 PSF using all the cameras. Recall that 
the DL07 model assumes that most of the dust is exposed to a single starlight intensity f/ min . 
A multi-pixel model with Um^j varying from pixel to pixel cannot be exactly reproduced by 
a model that assumes a single value of t/ min for the entire galaxy, hence fitting the summed 
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photometry with a single-pixel model will in general give a dust mass that will differ from 
that obtained by summing over a multipixel model. 

Figure [US shows that if all (IRAC, MIPS, PACS, SPIRE) data are used, the dust mass 
obtained by fitting the global SED is very close to that obtained from a multipixel fit of the 
resolved galaxy: the global dust mass is within 18% of the best estimate for NGC628, and 
within 3% of the best estimate for NGC 6946. If data from some cameras are not used, the 
single-pixel estimate for the dust mass will of course change, but we see that for NGC 628 
and NGC 6946 the resulting dust masses vary suprisingly little provided we have SPIRE data 
or MIPS160. If we use only PACS for A > 70yum, the dust mass estimate for NGC 628 is 
high by ~10%. The M160_lll_000_000 column in Figure [TBI compares the global dust mass 
estimated from a single-pixel fit using only Spitzer (IRAC and MIPS) data versus our "gold 
standard" best estimate (multipixel model using all cameras). The dust masses estimated 
using only the global photometry from Spitzer are within ~30 % of our present best estimates. 
Thus, at least for normal-met allicity star-forming galaxies such as NGC 628 and NGC 6946, 
dust mass estimates ba sed on IRAC and M IPS photometry, e.g., the dust masses estimated 
for the SINGS sample (IDraine et al.l 120071 ) appear to be reliable. 
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PSF and cameras used 

Fig. 17. — Dust mass estimates based on global photometry for NGC 628 and NGC 6946, using different 
camera combinations, relative to the sum of the dust mass obtained from pixel-by-pixel modeling of the 
resolved galaxy using the same camera set. 



Figure [T71 compares the global dust mass estimated using global fluxes to the dust masses 
estimated summing the dust masses over each map, for different sets of cameras. Essentially 
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it is a measure of the non- linear behavior of the dust modeling. For each camera and PSF 
combination, performing resolved modeling allows the dust to be colder in some regions, and 
this introduces more dust mass. Modeling based on global photometry under-predicts the 
dust masses by 0-20% in most cases. 

In principle, using the MIPS160 PSF results as our "gold standard" best estimate for 
the dust parameters may also suffer from the above-mentioned non-linear behavior: each 
MIPS 160 pixel covers a rather large area of the galaxy, in which the dust properties (e.g., 
<?pah) ot radiation field (e.g., ?7 min ) may have variations. Unfortunately, the discrepancy 
of PACS160 and MIPS160 fluxes makes it important to retain MIPS160 photometry, and, 
therefore, to use the broader MIPS160 PSF. 

Table H] summarizes the main parameters obtained by single-pixel and multi-pixel mod- 
eling of the galaxies at MIPS 160 resolution, using all the cameras. 

Table 4: Multi-pixel and single-pixel best-fit model parameters at MIPS160 resolution. 



Parameter 


NGC 628 
Resolved Global 6 


NGC 6946 
Resolved" Global 6 


M d (M ) 


2.87xl0 7 


2.33xl0 7 


6.74xl0 7 


6.62xl0 7 


L d (L Q ) 


6.83xl0 9 


7.08xl0 9 


3.31xl0 10 


3.20xl0 10 




1.46 


2.00 


3.05 


3.00 


u 


1.75 


2.23 


3.61 


3.55 


100 x / PDR 


11.6 


10.2 


14.2 


15.3 


100 x g PAH 


3.67 


3.70 


3.55 


3.30 


100 x M d /M H C 


0.82 


0.66 


0.63 


0.61 



a Total or mean values are defined in equations (1211 - 1241) . 
b Values are from a model fit to the global photometry. 
c For X co ,2o = 4. 



9. Alternative parameterizations of the starlight intensity distribution. 



Recently iGalliano et al.l (120111 ). in a resolved study of the dust in the Large Magellanic 
Cloud, advocated a starlight heating intensity distribution function different from what is 
e mployed in the c urren t work. They recommended using the starlight distribution function 
of iDale fc Heloul (120021 ). with a power law distribution between two adjustable parameters 
U min and U max , with adjustable power law exponent a. We note that IGalliano et al.l (1201 if ) 
uses a slightly different nomenclature: t/ max = U min + AU. The distribution used in the 
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present work (eq. [9] and [TO]) incorporates an additional "delta function" component at U n 
but fixes the value of C/ max = 10 7 . 



Galliano et al.l (120111 ) claimed that satisfactory fits can be obtained without the "delta 
function" componen t. We note that while we fix U max to a given value in our model, 
Galliano et al.l (1201 if ) u se Umax &s & fre e param eter. Both models have thus the same number 



of degrees of freedom. iGalliano et al.l (120111 ) based their claim on maps of the LMC using 
IRAC, MIPS and SPIRE. Here we test this claim using our maps of NGC 628 and NGC 6946. 

Figure [TBI shows maps of y 2 (goodness of fit) obtained for three different starlight distri- 
bution schemes, all done at MIPS160 resolution, using all 13 cameras. The top row images 
correspond to NGC 0628 and the bottom row images correspond to NGC 6946. In the left 
column the starlight distribution is a power law between U min and U max with power law 
exponent a (a ll adju sted) without the delta function contribution at U m - m , as proposed by 
Galliano et al.l ( 1201 ll ). It has 6 adjustable parameters (fi*, Md us t, <?pah, C^min, ^max, a), and 
therefore 13 — 6 = 7 degrees of freedom per pixel. The fit for this starlight distribution is 
poor, giving \ 2 m the range 10 — 20 in the bright regions of the galaxies. In the middle 
column the starlight distribution is a power law between U min and fixed U max = 10 7 , with 
power law exponent a (only U m i n and a are adjusted), plus a delta function contribution 
at f/min • This fit produces excellent results, and is the one recommended in the present 
work. We note that this fit has the same number of degrees of freedom as the ones in the left 
column panels; the adjustable parameter £/ max is replaced by the normalization of the "delta 
function" via the parameter 7. In the right column the starlight distribution is a power 
law between U m \ n and £/ max with power law exponent a (all adjusted), plus a delta function 
contribution at U m i n . Although the right panels have one extra adjustable parameter (£/ max ) 
with respect to the recommended fit, the fit quality is similar to the middle panels. If U max 
were a relevant parameter, one would expect \ 2 to decrease significantly (of order 1). In fact, 
in the right column fit, for most of the bright pixels the best-fit U max value is U max = 10 7 
and thus the value of y 2 is for those pixels exactly the same as in the center column fit. 

We found that the \ 2 values of the d istribution used in the current work are significantly 
smaller than the ones found using the IGalliano et al.l (1201 ll ) distribution. Therefore, we 
recommend adding the "delta function" component at U m i n . Since the fit does not improve 
significantly by allowing U max to be fitted, we recomme nd fixing it to 11™* = 10 7 , leading to 
a simpler, more robust fit. It should be noted that the IGalliano et al.l (120 111 ) study is based 
on the dust in the Large Magellanic Cloud, and in our present work we use NGC 628 and 
NGC 6946, two galaxies with metallicities closer to the Milky Way metallicity. Additionally, 
due to the proximity of the LMC, each pixel covers a physical area much smaller than our 
modeling pixels. 
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Fig. 18. — Comparison of \ 2 (goodness of fit) of three differ ent starlight d i stribu tion schemes. In the 
lef t column the starlight distribution is of the form proposed by iDale fc Heloul f|20Q2l ) and recently favored 
by iGalliano et al.l f|201 lh : a power law between U m [ n and Umax (both adjusted) without the delta function 
contribution at U m in- The fit for this starlight distribution is poor. In the center column the starlight 
distribution is of the form used in the present modeling: a power law between U m i n and U max = 10 7 , plus 
a delta function contribution at U m 'm- This fit produces excellent results. In the right column the starlight 
distribution is a power law between U m [ n and /7 max (both adjusted), including a delta function contribution 
at U min . Although the right panels have one extra degree of freedom (/7 max ), the fit quality is similar to the 
middle panels. The top row is for NGC 628 and the bottom row is for NGC 6946. The power law exponent 
a is treated as an adjustable parameter in all the fits. 
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10. Discussion 

10.1. Dust Mass Estimation 

Dust mass estimates are, of course, model-dependent. Above we have estimated the dust 
mass using a specific dust model, the DL07 silicate-graphite-PAH model, with an assumed 
parametric form for the starlight intensities heating the dust. The physical dust model 
and the ansatz for the distribution of starlight intensities together appear to successfully 
reproduce the observed SEDs, and to give reliable estimates of the dust masses. 

As discussed in §4j the far-infrared opacity of the amorphous sil icate originally put for- 



ward by DL84 was subsequently adjusted slightly ( iLi fc Draindl200ll ) to improve agreement 
with the far- infrared and submm emission observed for the local high-latitude dust. Thus 
the DL07 model uses dust properties that can reproduce the observed FIR-submm emission 
from the Milky Way (MW) cirrus with a single starlight heating intensity - no "cold dust" is 
needed for the MW cirrus. In this dust model, the opacities are assumed to be independent 
of the dust temperatures. 

Here we see that the same dust model, with suitable adjustment of the starlight assumed 
to be heating the dust, is able to reproduce the emission from NGC628 and NGC6946 out 
to 500/im, without introducing any "cold dust" component. The modeling performed at 
MIPS160 resolution, using all the IRAC, MIPS, PACS, and SPIRE cameras, in addition 
to being able to reproduce the observed SED, gives dust masses that are in line with the 
expected dust/H mass ratios for these galaxies: the dust/H mass ratio images in Figures E] 
and[T0l at MIPS160 resolution, are smooth, and the global dust/H mass ratios are 0.0082 ± 
0.0017 and 0.0063 ± 0.0009 for NGC 628 and NGC 6946, respectively. 

Observed de pletions in the local Milky Way indicate a dust/H mass ratio M&^/M-r = 



0.0091 ± 0.0006 ( lDraindl2011bl . Table 23.1). Thus, a galaxy with a similar fraction of inter- 



stellar heavy elements in dust would be expected to have 



M, 



dust 



M H 



0.0091 x 10 [o] , (29) 



where [O] = log 10 [(O/H)/(O/H) Q ]. Thus galaxies with heavy element abundances (and 
depletions) similar to the local Milky Way should have M^t/Mn ps 0.009. 

NGC 628 and NGC 6946 appear to be mature star-forming galaxies that would be 
expected to have interstellar metallicities similar to the local Milky Way. If so, then the 
values of Md us t/^H found by fitting a dust model to the observations appear to be in excellent 
agreement with e xpect ations. We note that the HII region oxygen abund ances found by 



Moustakas et all ( j2010h together with (A Q ) Q = 12 + log lo (O/H) = 8.73 flAsplund et al. 
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2009|) yield [O] = -0.46 and -0.36 for NGC 628 and NGC 6946, and Eq. ([29]) would predict 
M d /M H = 0.0032 and 0.0040 for NGC 628 and NGC 6946, respectiv ely. However, we s uspect 
that the PT05 oxygen abundance estimates may be biased low: iMoustakas et al.l (120 lOf ) 
list PT05 HH-region oxygen abundances for 38 galaxies; the highest oxygen abundance is 
Ao = 8.59 ± 0.11 for NGC 4826. It seems unlikely that none of the galaxies in their sample 
have oxygen abundances that are solar or supersolar. 



10.2. Dust Opacity in Molecular Gas and the Value of X, 



CO 



Both NGC 628 and NGC 6946 are rich in molecular gas, and the estimated gas mass 
depends on the adopted value of Xco- In the present study we have adopted Xco, 20 = 4, 

2 found 



which, as discussed in Section 12. 3[ is significantly larger than the value Xco, 20 
for resolved CO clouds in the Milky Way. The larger value of Xco adopted here may reflect 
the presence of so-called "dark gas", diffuse H 2 with very l ow CO abundances , which does 
not r adiate effectively in either H I 21-cm or CO J = 2 — > 1 (jWolfire et al.ll2010t iLeroy et al. 
201lh . 



However, if the dust opacity in molecular regions is actually larger than the opacity in 
H I gas, then the present approach - where we favor a value of Xco that minimizes small- 
scale structure in maps of dust op tical depth/H surface density - will overestimate X co . 
Planck Collaboration et al.l (j2011al ). usin g measurements of submm emission by Planck, and 
iVn inferred from NIR reddening of stars (IPineda et al.ll2010l ). conclude that the dust opacity 
per H nucleus in the Taurus molecular cloud is larger than in the local diffuse H I by a factor 
R m 2.0 ± 0.4. Such an enhancement in the far-infrared and submm opacity might be a 
consequence of coagulation, which is expected to incr ease the far-infrared and submm opacity 
( lOssenkopf fc Henninglll994l ; IStognienko et al.lll995l ). However, grain coagulation could also 
flatten the NIR extinction curve, so that the value of Nn inferred from the (J — H) and 
(H — K) stellar colors might be an underestimate. It should also be noted that studies of 
the Corona Austrina molecular cloud, using MIPS160/LABOCA870/im ratios to determine 
the dust temperat ure, found a normal ratio of 870/im optical depth to visual extinction 
iluvela et al]bo09h . 



The dust model used here has been calibrated on dust in H I regions. If the actual dust 
opacity per H nucleon in molecular clouds is larger than in H I clouds by a factor R, then the 
actual value of X co in NGC 628 and NGC 6946 would be X co ,2o = ^/R- A value of R « 2 
would then bring us into agreement with the Xco, 20 ~ 2 inferred from other estimators 
(virial mass estimates, 7-ray emission) of molecular mass. 
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10.3. Dust Mass Estimates from Single- Temperature Fits 



Skibba et al.l (120111 ) used a single dust temperature T d with an assumed k v oc A 1,5 
opacity to fit the 70-500/im photometry from MIPS and SPIRE. The opacity at 500/im was 
taken to be that of the DL07 dust model. They estimated T d = 24.0 K, M d = 1O 7 O3±O O8 M 
for NGC 628, and T H = 26.0 K, M d = lO 7 - 47±ao8 M for NGC6946. The dust masses esti- 



mated by ISkibba et al.l ( 1201 ll ) for these two galaxies are smaller than the dust masses found 
here (see Table |3]) by factors of 3.3 and 3.6, respectively. This is the result of using a single 
dust temperature to try to reproduce emission from 70-500/im. With the more realistic 
assumption of a distribution of dust temperatures, a small amount of warmer dust can pro- 
vide much of the 70/im emission, thus requiring an increased mass of "normal" temperature 
dust to account for the emission at A > 160//m. A range of dust temperatures is of course 
expected from both spatial variations in the starlight intensity heating the dust, and the fact 
that a grain model with more than one dust composition, and a broad range of grain sizes, 
will have a distribution of temperatures even in a single radiation field. Single-temperature 
dust fits, if constrained by emission at 70//m , will not provide a reliable estimate of the dust 
mass. In a separate work (Aniano & Draine 2012, in prep.) we discuss the bias introduced 
when the DL07 SED is approximated by a (single or dual) temperature blackbody multiplied 
by a power law opacity. 



10.4. Radial Gradients in Dust and Starlight 

In both NGC 628 and NGC 6946, we find that the dust/H mass ratios outside the 
nucleus vary slowly, with little indication of a decrease in the dust/H ratio as one moves 
outward (see Figures HI and|9j). In NGC 6946, however, the dust/H mass ratio appears to 
have a pronounced minimum at the center. We interpret this as due to overestimation of 
the gas mass in this r egion: we employ a single value of X co ,2o — 4 for this galaxy, but 



Donovan Meyer et al.l ( 120121 ). using virial mass estimates, find Xqq^q = 1.2 for the giant 
molecular clouds (GMCs) in the central 5 kpc. Because the molecular gas dominates near 
the center, our use of a higher value of Xqo implies an overestimate of the gas mass, resulting 
in an underestimate of the dust/H ratio. We therefore suspect that the central minimum 
in the dust/H ratio in Figure [TD]is entirely an artifact of using a value of Xqq that is too 
large for the central region. If Xco,2o = 1-2 is appropriate in the center of NGC 6946, 
the gas surface density will be lowered by a factor ~ 1-2/4, and the dust/H mass ratio 
increased by a factor ~ 3, from the central value ~ 0.004 in Figure 7 to the expected value 
~ 0.012. The present study of the dust surface density therefore supports the finding by 



Donovan Meyer et al.l (120121 ) of a low value of Xqq in the center of NGC 6946. 
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In both NGC 628 and NGC 6946, the PAH abundance parameter (/pah appears to be 
quite uniform, with little evidence for a radial decline (see Figuresl5land [T0l). Previous studies 



have shown that low-m etallicity galaxies have low values of (/pah (e.g., lEngelbracht et al. 
20051 ; iDraine et al.ll2007l ); a ga laxy with a negative radial g radient in metallicity might then 
show a radial decline in (/pah- iMunoz-Mateos et al.l (120091 ) found radial declines in (/pah for 
a number of SINGS galaxies. For both NGC 628 and NGC 6946, (/pah appeared to be quite 
uniform out to a radius of ~10kpc, in qualitative agreement with our findings here. Future 
papers will examine radial trends in the KINGFISH sample. 

The starlight intensity parameter ?7 m i n is presumed to characterize the average starlight 
intensity in the diffuse interstellar medium. The maps of U m i n (Figs. |5ghi and ITOghi) have a 
peak U-mir, ~ 4 and 6 in the central regions of NGC 628 and NGC 6946, respectively, declining 
to U m i n ~ 0.3 near the edge of the galaxy mask (~10kpc from the center). It is gratifying 
that ~ 1 at ~8 kpc from the center, just as for our location in the Galaxy. 

The parameter /pdr is the fraction of the dust heating that is produced by starlight 
intensities U > 100. Averaged over the full galaxy, (/pdr) = 0.12 and 0.14 for NGC 628 and 
NGC 6946, respectively, but both galaxies have hot spots where /pdr is much higher, with 
peak values of ~0.30 (see Figs. EJkl and[l0]kl). With the ~635 pc resolution of the SPIRE250 
camera at the 7.2 Mpc distance of NGC 628, these hot spots presumably correspond to 
regions of very active star formation. 



10.5. Interpretation of the Starlight Heating Parameters a and 7 

As discussed in §4.4} a fraction 7 of the dust mass is taken to be heated by a power-law 
distribution of starlight intensities between U m i n and U max with dM/dU oc U~ a , and the 
remaining fraction (1 —7) of the dust mass is heated by a radiation field with intensity U m i n . 
Dust heated by this simple parameterization of the starlight intensities is quite successful in 
reproducing the observed SED (see the SEDs in Figures [7b,f; [Sb,f,i,l; IT2"b.f; and[T3"b.f.i.l). 



Dale et al.l (120011 ) discussed two ideal cases that have analytical predictions for the value 
of a: a single star in a homogeneous diffuse medium, and a dark cloud. The first case, with 
U oc r~ 2 , has dM d /dU oc U~ 2 ' 5 , i.e., a = 2.5. In this case, U max would be very large, > 10 10 , 
the heating rate at which dust grains would vaporize. In the second case, a slab where the 
heating intensity is primarily attenuated by dust absorption, one would get dM d /dU oc U -1 , 
i.e., a — 1, and U maK would be set by the starlight intensity at the edge of the slab. This case 
might correspond to a photodissociation region, the edge of a dark cloud. Different clouds in 
the pixel would have different values of U max , so their superposition could be approximated 
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Fig. 19. — Characterization of the starlight power-law component, at the resolution of PACS160 (left), 
SPIRE250 (center), and MIPS160 (right). Row 1 and 3: power-law exponent a map for NGC 628 and 
NGC 6946 respectively. Row 2 and 4: dust mass fraction in the power-law component 7 map for NGC 628 
and NGC 6946 respectively. 
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as a power-law with a slightly larger value of a. We expect dust in a real galaxy to have 
1<«<2.5. 

Consider for the moment the fraction 7 of dust with dMj dU oc U~ a . The power radiated 
is dL oc U dM oc U 2 ~ a d log U. A value of a = 2 would have uniform dust luminosity per 
unit interval in log U. A value a < 2 would concentrate the dust luminosity L in the high 
radiation field regions (U ~ £/ max ), and a value a > 2 would have L dominated by dust with 

Figure [191 shows the starlight power-law component parameters a and 7. The power-law 
index a has a preferred value a ~ 1.6 in the bright regions. In these regions, the power-law 
distribution is representing the heating of dust in high-f/ regions, e.g., photodissociation 
regions near OB stars. There are also regions with a ~ 2.4. In these regions, most of the 
dust luminosity is concentrated in regions with U ~ U m i n , and the power law component is 
essentially making the U = Umi n component broader, i.e., is allowing for variations in the 
starlight intensities in the diffuse ISM. 

We do not attach great physical significance to the parameters 7 and a. The more 
physically meaningful quantities are the mass-weighted mean starlight intensity U, and /pdr, 
the fraction of the dust luminosity that originates in regions with U > 10 2 . In Figures [5] and 
[TUl one sees that /pdr > 0.05 over most of the galaxy mask for both galaxies. 



11. Summary 

Using Spitzer and Herschel data, we perform a resolved study of the dust physical 
parameters, and of the starlight heating the dust, in the nearby galaxies NGC 628 and 
NGC 6946. We employ the DL07 dust model, having amorphous silicate grains and car- 
bonaceous grains, including PAHs (see $4]). The model has a distribution of grain sizes, 
and we allow for a distribution of starlight intensities heating the dust within each pixel. 
The model assumes a frequency-dependent opacity that reproduces the observed emission 
spectrum of high- latitude dust. 

In order to perform resolved studies of the galaxies, it is important to convolve all the 
images into a comm on PSF. This is achieved using the convolution kernels generated by 



Aniano et al.l (120111 ). We perform our dust modeling using all appropriate combinations of 



PSFs and cameras. Table [T] lists some of the resolutions and cameras used in this work. 
Our principal findings are as follows: 



1. The dust model is quite successful in reproducing the observed SED over the full 
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wavelength range 6-500/iin where dust emission dominates (see the SEDs in Figs. [7b, f; 
I8"b.f.i.l; [T2b.f; and [T3b.f.i.l; and the model/observed maps in Figs. [6b,f and ITTb.f). The 
dust model reproduces the emission out to 500yum without introduction of a "cold dust" 
component, and with no allowance for possible temperature-dependent dust opacities. 
The DL07 dust opacities therefore appear to be consistent with both the observed 
emission from local "cirrus" at high galactic latitudes, and the observed SED from 
500pc-sized regions of the ISM in NGC 628 and 6946. There is no indication that 
T-dependent opacities are needed to reproduce these data. 

2. Maps of the dust/H mass ratio show it to be relatively uniform over both galaxies, 
outside of the nucleus (see Figs. HI and0). The derived dust/H mass ratio for these 
galaxies is consistent with that expected if the interstellar abundances in NGC 628 
and NGC 6946 are close to solar. As near-solar abundances seem likely, this is strong 
support for the quantitative accuracy of the dust masses obtained by modeling the IR 
and FIR emission. 

3. Figure [141 shows how dust mass estimates depend on the camera and PSF used. The 
"gold standard" is taken to b e dust modelin g done at MIPS 160 resolution, using all 



cameras, with Scanamorphos flRousselll2012l ) processing of the PACS data. Relative 



to this standard, mass estimates done with smaller pixels (giving up the lowest reso- 
lution cameras) tend to be biased slightly high. At SPIRE250 resolution, the bias is 
~38%. Working at SPIRE350 resolution provides a good compromise between resolu- 
tion and accuracy, with ~30% errors for the total dust mass estimated for NGC 628 
and NGC 6946. 

4. Even though the dust modeling process is non-linear, resolved modeling is compatible 
with modeling using the global photometry of the galaxy. Differences in the inferred 
parameters are in most cases under 20%. 

5. In NGC 6946 the dust/H mass ratio calculated with a single value of ATco,2o = 4 appears 
to have a minimum at the center (see Fig. ITOk.b.c). This minimum is interpreted as an 
artifact due to overestimation of the molecular mass. The dust surface densities found 



here therefore support the finding by iDonovan Meyer et al.l (120 121 ) of Xqo,2o ~ 1-2 in 
the center of NGC 6946. 

6. The present fitting procedures employ six adjustable parameters for each pixel: Q±, 
M dust , gpAH, U min , a, and 7 (we fix U max = 10 7 ). For NGC 628 and NGC 6946 we 
fin d that this approach p rovides substantially better fits than the approach advocated 



P 





bv iGalliano et al.l (120111 ) who treat Q+, M dust , (/pah? ^m in, a, and U max , as adjustable 
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parameters and fix 7 = 1: compare the \ 2 maps in Fig. [T8k with H~8b . and [T8H with 

EEfe. 

7. The starlight heating the dust within a single pixel is characterized by three adjustable 
parameters: U m i a , 7, and a (see eq. |9]and[T0]). U m i n is interpreted as the intensity of 
the diffuse starlight, responsible for the bulk of the dust heating. Maps of U m i n (see 
Figs. 2 and 7) show significant structure: there is a tendency of U m i n to be higher in 
spiral arms, as well as a significant radial decline in U m \ n . The overall mass-weighted 
mean value of U m i n = 1.5 for NGC 628 and U m i n = 3.1 for NGC 6946, but U m i n declines 
to values of 0.3 or lower in the outer regions of both galaxies. 

8. The parameter /pdr is the fraction of the dust heating power that is contributed by 
starlight intensities U > 100. The overall value is (/pdr) = 0.12 for NGC 628 and 
(/pdr} — 0.14 for NGC 6946, but in both galaxies /pdr peaks at local maxima of the 
dust surface brightness, which are associated with star-forming regions. At SPIRE250 
resolution, /pdr reaches values as high as 0.25 in both galaxies. 

9. Current PACS photometry disagrees with MIPS photometry (see Fig. [27] and 1281) . 
PACS70 photometry is up to 80% higher than MIPS70, and PACS 160 is up to 50% 
higher than MIPS70 photometry in the bright areas of the galaxies. 

10. The PACS-MIPS photometry disagreement induces a bias when trying to model dust 
at high spatial resolutions (PACS160 PSF), where MIPS70 and MIPS160 cannot be 
used. We do not recommend modeling dust at PACS 160 resolution in low surface 
brightness areas. Modeling done at SPIRE250 resolution is more reliable than modeling 
at PACS160, but the inferred dust masses still disagree with our "gold standard" (i.e., 
modeling done at MIPS160 resolution using all the IRAC, MIPS, PACS and SPIRE 
cameras) by up to ~ 30%. 



Subsequent work (Aniano et al. 2 012, in preparation) w ill extend this study to all 61 
galaxies in the full KINGFISH sample f kennicutt et al.lboill ). 

We are grateful to R.H. Lupton for availability of the SM graphics program, and to the 
anonymous referee for helpful suggestions. 

This research was supported in part by JPL grants 1329088 and 1373687, and by NSF 
grant AST1008570. 
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A. Separation of the Target Galaxy and "Non-Background" Regions from 

Background Regions 

In order to backgound-subtract an image, one needs to identify the image areas with- 
out any emission, i.e., recognize all the "non-background" areas to exclude them in any 
background estimation procedure. Different cameras have different sensitivities to extended 
emission and other sources, so if one attempts to background-subtract an image using only 
the image itself, one may under- or over-overestimate the background level. As an example, 
the PACS cameras are not very sensitive to the extended emission in the North part of 
NGC 6946, (see Fig. [2]), so, unless these areas are recognized by other cameras, one could er- 
roneously consider them as background areas leading to over-estimation of background level. 
We proceed using all the information available: i.e., using all the cameras, to recognize the 
"non-background" areas 1^1 

We developed an algorithm (described in Appendix B) that, given an image / and a 
set of "non-background" areas, identifies "background" areas and estimates the variance in 
the "background" area (i.e., a measure of the image noise). To obtain the most accurate 
background estimation, we generate background masks for each camera and combine them 
to generate a mask of the non-background areas, proceeding in several steps as follows. 

1. As a preliminary identification of the "non-background" areas we select the sky regions 
with S/N > 2 in the most sensitive cameras. In our study, for this purpose we use 
IRAC5.8, IRAC8.0, MIPS24, MIPS70, and MIPS160. We estimate the noise in the 
images using the iterative algorithm described in Appendix B, using an empy set as 
starting "non-background" areas. The "non-background" regions obtained contain 
the target galaxy and bright background sources. We consider the Starting Non- 
Background Mask (SNBM) as those regions in which at least 3 of the 5 cameras have 
S/N > 2. The SNBM include the target galaxy and some foreground bright sources. 

2. For each band b, we obtain a Preliminary Individual Background Mask (PIBMfc) with 
the iterative algorithm described in Appendix B. For this construction, we consider 
the pixels in the SNBM as "non-background" pixels (i.e., they were not considered as 
background pixel candidates). 

3. We construct a Preliminary Background Mask (PBM) as follows: We compute the 
average A of the PIBM&s, where only the cameras observing the pixel are taken into 



1 Non-background" areas will include the target galaxy, and any other recognizable emission structures 
in the field 
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account. At each pixel, the value of A will be the fraction of cameras that do not detect 
any emission with S/N > 1.8 x 0.4202 = 0.76 . The outer regions of the field (that may 
be observed by a small number of cameras), can still be considered as background if the 
observing cameras do not show emission. We smooth the image A by convolving it with 
a (normalized) gaussian kernel oc exp[— (r/ro) 2 ], with tq equal to 0.5 times the width of 
the pixel. The convolution lowers the value of A near sources, reducing the background 
areas. We consider the PBM as the pixels where the smoothed A > 0.35, i.e., fewer 
than 35% of the cameras detect emission. With 13 cameras, PBM corresponds to the 
areas with S/N > 0.76 detection in no more than 4 cameras. 

4. Using the complement of PBM as "non-background areas" and the algorithm described 
in Appendix B, we background-subtract each individual image. The complement of 
PBM masks all recognized sources in the background estimate for each individual 
image. 

5. Using the background-subtracted images, we proceed to perform a preliminary dust 
modeling at MIPS160 PSF (FWHM = 38.8") using all the cameras available. We 
obtain a preliminary estimate of the dust luminosity surface density S ^ d in each pixel. 
We will construct a Preliminary Galaxy Mask (PGM) based on the dust luminosity 
surface density obtained. In order to obtain smoother galaxy borders, we smooth 
the E^ d image by convolving it with a (normalized) circular kernel of 20" radius. We 
consider the PGM as the regions with (smoothed) Ex, d > Si dimin . We use the threshold 
values £ Ld , min = 3.14 x 10 6 L o kpc~ 2 for NGC 628 and £ Ld ,mm = 1.12 x lO 7 L kpc~ 2 
for NGC 6946. The thresholds S^ d min were manually chosen so that the FGM covers 
the area where can reliably estimate gpAH, U min , 7, and a, i.e., so that the inferred 
dust parameters have S/N > 1. Regions not connected to the main galaxy (i.e., other 
sources) are now omitted from the mask. Figure |2U] shows the histogram of the dust 
luminosity Ei d for a region containing NGC 628 (left) and NGC 6946 (right). The 
solid shaded (red) area corresponds to the pixels identified as "galaxy pixels" (i.e., in 
the galaxy mask), and the line-shaded (blue) pixels correspond to "non-galaxy pixels" 
(i.e., background pixels and other sources in the field not connected to the galaxies). 
We note that the galaxy cut-off is performed in the smoothed image, and thus, the 
boundaries in the non-smoothed luminosity histograms shown in Figure |20] are not 
sharp. 



At this point, we have a well-determined separation of the image into (preliminary) 
"galaxy pixels" and the non-galaxy pixels. We now continue with steps 6-9, where we 
essentially repeat steps 2-5, starting with the PGM as our SNBM. The new masks generated 
in steps 6-9 are very similar to the corresponding ones previously generated on steps 2-5, 
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Fig. 20. — Dust luminosity surface density Si d histogram of a region containing NGC 628 (left) and 
NGC 6946 (right) shown in Figs. 4. and 10. The modeling is done at MIPS160 resolution (18"pixcls) using 
all the cameras available. The shaded (red) area corresponds to the pixels identified as "galaxy pixels" in 
the galaxy mask. The line-shaded (blue) area corresponds to "non-galaxy pixels" (i.e. background pixels 
and other sources in the field not connected to the galaxies). 



with only small differences in the low surface brightness areas. The combined procedure is 
robust and precise. 

6. For each band b, we obtain an Individual Background Mask (IBM&) with the iterative 
algorithm described in Appendix B. For this construction, we consider the pixels in the 
PGM as "non-background" pixels (i.e., they were not considered as background pixel 
candidates) . 

7. We construct a Background Mask, (BM) using the same procedure as in step 3, but 
using the IBM b s instead of the PIBM fc s. 

8. Using the the complement of the BM "non-background areas" and the algorithm de- 
scribed in Appendix B, we background-subtract each individual image. The comple- 
ment of BM masks all recognized sources in the background estimate for each individual 
image. 

9. Using the background-subtracted images, we proceed to perform the final dust mod- 
eling at MIPS160 PSF using all the cameras available. We construct a Galaxy Mask 
(GM) in the same way as we constructed the PGM (using the same values of T<L d ,min)- 
The final GM and PGM will differ only in a few boundary pixels. 
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In Aniano et al. (2012, in preparation) the same masking method is applied to the full 
sample of 61 KINGFISH galaxies. Background companion galaxies are manually identified 
and removed from the masks when possible. 

B. Algorithm for Background Recognition (Image Segmentation) and 

Subtraction. 

Given an image I expressed in its native pixel grid (x,y), and a (potentially empty) 
set NBstd of "non-background" sky regions expressed in a standard 4" x 4" grid (a grid 
common to all the images), we iteratively construct a sequence of "background" pixels masks 
B n ,n = 0,1,2,..., and a sequence of tilted planes plane n ,n = 0,1,2,... that approximate 
the image values over the sets B n ,n = 0,1,2,.... The sets B n ,n = 0,1,2,... will exclude 
the "non-background" NBstd regions, i.e., it is assumed that the pixel in NBstd are not 
background-candidate pixels. 

Throughout the algorithm, the given "non-background" regions will be masked and not 
considered in the fitting. We start by bi-linearly interpolating or averaging (depending on 
the grid relative pixel sizes) the set NBstd into the image / grid (x,y) to obtain NBp rev . 

We first smooth the original image by convolving it with a (normalized) gaussian kernel 
oc exp[— (r/r ) 2 ], with r equal to 0.56 times the width of the native pixel for each camera. 
The factor 0.56 is chosen so that each pixel contributes half of its flux to the corresponding 
pixel in the convolved image. This step provides a less noisy image, in which we can more 
robustly identify the galaxy and the external sources. 

We proceed as follows. We start by setting B to be all pixels in the grid (x, y), except 
for the pixels in NB prev (if any). We fit a plane planeo to the image values over the -Bo pixels. 
This is achieved by finding the planeo that minimizes eo defined as: 

e ™ = [I(x,y) - plane n (x,y)} 2 , (Bl) 

where n — 0. We compute the dispersion #0 defined as: 

° n s vW,y)e£n}-^' (B2) 

where n = 0, so the sum is over the B (original "background" candidate) pixels. 

We take the first (n = 1) set of "background pixels''^ B\ to be all pixels that do not 



"Background pixels" consist of pixels that are not part of the target galaxy, and any other recognizable 
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belong to NB prev , or lie more than k8q above planeo (k is a constant, and we use k = 1.8, 
see below for details): 



In the first step, these pixels typically correspond to the brightest point sources and central 
pixels of the galaxy. 

We now proceed iteratively by: 

1. Fitting a tilted plane plane n to all the candidate "background pixels" B n , by minimiz- 
ing e n defined by eq. fIBlj) . 

2. Computing n , defined by eq. flB2j) . where the sum extends over the B n pixels. 

3. Identifying a new set of "background pixels" B n+ i to be those pixels which lie inside 
the interval (plane n — n6 n ,plane n + K,8 n ), and do not belong to NB prev : 



This iteration will generate a sequence of tilted planes plane n , and a (generally decreas- 
ing) set of "background pixels" B n . The iteration will converge, i.e., it reaches a step M 
where B M+1 = B M ), typically in M = 20 — 250 iterations. 

Once the iteration converges, we have set of "background pixels" Bm, and a fitting plane 
plantM (all the computed sets B n and planes plane n will depend on the choice of n, see below 
for a discussion). We bi-linearly interpolate or average Bm to the standard 4" x 4" grid. This 
will generate a mask of the "background regions" in the standard grid, so it can be compared 
with the results of the algorithm from images of different cameras. We can generate standard 
mask of "non-background regions" by taking the complement of the "background regions" 
mask. The left column of Figures [22] and [23] show the resulting residuals for NGC 628 and 
NGC 6946 respectively over B M , for k = 1.8. Rows 1-4 correspond to IRAC4.5, MIPS24, 
PACS160, and SPIRE350, respectively. It can be seen that the algorithm masks out all the 
point sources, extended sources and several background regions. Residuals slightly positive 
are found in the edges of the extended emission areas, although most of the extended emission 
areas are recognized and excluded. 



emission structures (e.g., background galaxies) in the field, removing also any pixels that are "outliers", 
cither above or below the general background. 





(B4) 
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After the iteration converges to a set of "background pixels" Bm, if we desire to perform 
a background subtraction to the images, we run the iterative steps again. In the second 
iteration, we use the original (non-smoothed) image, and the complement of the set Bm 
found as the starting set of "non-background pixels" NB prev , i.e., the pixels which were 
recognized as having emission in the first iteration are masked out. The second run of the 
iteration will converge to a new set of "background pixels" Bn, and a new fitting plane 
plane n, typically in N = 20 — 250 iterations. Clearly, by construction, we always have 
Bn Q Bm- Finally, the last plane plane^ is subtracted from I. 

The reason to perform the second iteration, with the original image and a larger set 
of starting "non-backgrund pixels" is so we can remove noisy pixels from the background 
pixels, recognize the edges of the emission areas more precisely, and perform a more precise 
background level estimation. Moreover, if we try to perform a single iteration (i.e., by going 
directly to the second iteration) , then the larger image noise will prevent us from recognizing 
the dimmer sources and extended emission areas. The center column of Figures [22] and [231 
show the resulting residuals for NGC 628 and NGC 6946 respectively, over the sets Bn- 

It can be shown that if there exist an infinite set of background points B in which the 
intensities are independent noise drawn from a distribution with probability density p, then 
the algorithm will asymptotically converge: B n > Bn Q Bo, that will depend on re. We 

n— >oo 

define /(re), the fraction of the original background-candidate pixels that are asymptotically 
considered background: 

/(re) = #{B N }/#{B } (B5) 

For noise distribution function densities p which are analytical around their expectation 
values, it can be shown that if re < \/3 then B n will converge to an empty set 0; and B n will 
converge to a non-empty set provided that re > v3, i.e., /(re < a/3) = 0, /(re > a/3) > 0. 
In particular, for a top hat distribution (p(x) = 1/2 -<=>- \x\ ^ 1), B N = B if re > a/3, i-e., 
/(re > V3) = 1- If the noise has a gaussian distribution with standard deviation a, we can 
numerically compute /(re), and the actual a/9 (the ratio between the distribution standard 
deviation a and the computed final dispersion 9). Figure [21] shows the values of /(re) (dotted 
green line) and 9 jo (solid blue line) for 1.6 < re < 2.6 for a gaussian distribution. The vertical 
red lines correspond to the threshold value re = \/3> and the used value re = 1.8 

Our choice of re is dictated by (1) the need to have re > a/3 in order for the procedure 
to converge to a non-empty set, (2) the desire to minimize inclusion of outliers that may be 
due to rare events (e.g., cosmic rays, image artifacts), and (3) the desire to be as sensitive 
as possible to dim sources. We choose re = 1.8 as a good compromise. If the background 
noise were gaussian, for our choice re = 1.8, we have /(re = 1.8) = 0.5506, 9/a — 0.4202, 
(i.e., the final iteration should fit only m 55% of the background points, and their measured 
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dispersion should be 0.4203 a). 

Algorithm performance for Gaussian noise 

1 

0.9 
0.8 
0.7 
0.6 
0.5 
0.4 
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0.2 
0.1 
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K 

Fig. 21. — Theoretical performance of the iterative algorithm for gaussian random noise. The horizontal 
axis correspond to the cutoff threshold k. The dotted green line shows 0/cr, the ratio between the dispersion 
of the final set Bn and the dispersion of the original distribution. The solid blue line shows the fraction of 
the original points that survive in Bn. The vertical red lines correspond to the critical value k = \/3 and 
the used value n = 1.8. 

How well does the iterative procedure behave with real images? Typically the iterative 
procedure converges in 20-250 steps. The right column of Figures 1221 and [231 show the result- 
ing histograms of the background-subtracted intensities in Bn for NGC 628 and NGC 6946 
respectively. The solid red curve is a gaussian distribution function with o = 0/0.4202, i.e., 
with the standard deviation expected from the measured dispersion. The fit of the gaussian 
curve to the histograms in Figures [22] and [23] is very good, supporting the validity of the 
presented algorithm. 

This background recognition algorithm has several advantages over other approaches. 
Firstly, the algorith is very robust, automatic, and quick, and produces superior results when 
dealing with all the ~ 1400 images of the KINGFISH sample. Secondly, by keeping only a 
subset (~ 55%) of the background pixels, the last computed dispersion is small enough that 
even dim sources are easily recognized and masked. This allows us to exclude virtually all 
detectable foreground or background sources from the subset of points used to determine the 
background of the images. Thirdly this algorithm also recognizes areas with dim extended 
emission. Fourthly, image artifacts over the background areas are easily identified, allowing 
us to minimize bias induced by them in the background estimation. This allows us to avoid 
over- or under-subtraction of the background in each image, which is critical in order to 




model the galaxy low surface brightness areas. 



Table 5: Background information for the different Cameras 
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Number of native pixels that are considered "background" . 



Percent of the non-galaxy pixels that are considered "background" pixels. 
c Standard deviation of the Gaussian fit of the "background" pixels. 



C. Correction for Missing Data (Bad Pixels) in the Original Images 

Some of the original images include bad pixels due to detector saturation, incomplete 
scanning, or other problems. This effect is not important when present in the pixels far away 
from the target galaxy, but may have undesired consequences if not treated appropriately 
when present near or within the target galaxy. A similar situation arises when observations 
are made over a small part of a galaxy, for example with the IRS instrument onboard the 
Spitzer Space Telescope. We develop a method for correcting this problem. 

In order to estimate (and correct) the influence of bad pixels of a camera A, we start 
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Fig. 22. — Analysis of the residuals in NGC 0628. Row 1: IRAC4.5. Row 2: MIPS24. Row 3: PACS160. 
Row 4: SPIRE350. Left column: residuals on the set Bm (after the first iteration of the algorithm on the 
smooth image). Center column: residuals on the set Bn (after the second iteration on the original image). 
Right column: histograms of the pixels on Bm- The (red) solid line is a gaussian with a = 0/0.4202 (see 
text for details). 
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3 but for NGC 6946. The bright carbon star V0778 Cyg, located near 
RA=309.044, Dec= 60.082 (slightly off the bottom left edge of the maps shown) produces strong artifacts 
in the IRAC and MIPS cameras, reducing the useful background area of the images. 
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by choosing a second camera B. We require that the sky coverage of camera B be larger 
than camera A, or at least that does not present bad pixels in the same areas as A. From 
all the possible candidates for being the correcting camera B, we choose the camera with 
wavelength most similar to A. For example, we typically correct the MIPS 160 camera with 
the SPIRE250 camer£0. 

We divide the original image grid of the camera A, OGa, into two set of pixels: OGA, goo d, 
the pixels of OGa with good A data, and OGA,bad, the pixels of OGa m which the camera 
A has bad data. 

We interpolate the background-subtracted image B from its original grid OGb into the 
OGa grid, obtaining B interp A, which looks morphologically similar to A. We construct 2 
images B good and B bad , both defined on OGa, as follow: 

B good (x, y) G OGa — I ^ nterpA ^ ^ V \ I °^°° d (CI) 

I [X, y) G OG A ,bad 

r> I \ r- r,n J (X,y) E OG A ,good , no x 

B had {x,y) eOG A = < , ,^ nn (C2) 

L BinterpA "<=>• [X , y ) G U(j A,bad 

We process the images B good and exactly as we do with the image A; we rotate 
them, convolve with the kernel that corresponds to the image A, and interpolate them to the 
final grid (FG). A given output pixel in FG will have flux flux good B, and fluxbadB in the 
final stage of processing the B good and B bad images respectively We define Fa,b(x, y) G FG, 
the fraction of the total flux in the final grid coming from areas with good A data, as: 

flux good B rnn\ 

flUX good B + flUXbadB 

Fa,b provides an estimate of how much flux we are missing in each output pixel due to 
the missing data in the original image A. In the pixels for which Fa,b > 0.65, (i.e, those 
where most of the flux is actually coming from regions where we know A) we can further 
compensate the final stage of the A image processing by multiplying the value obtained by 
1/F. In the points Fa,b < 0.65 a significant amount of the flux (more than 35%) is coming 
from regions where A is not known, so we mask and ignore them. 



19 We avoid correcting the MIPS160 camera with the PACS160 camera due to the inability of PACS160 
camera to correctly measure the flux in extended emission regions. Moreover, the SPIRE maps extend to a 
larger area than the PACS maps. 
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This correction algorithm would be exact if the color A/B were constant, and the PSFs 
of both cameras coincided. We test the algorithm by purposely removing some part of an 
image, and correcting it with itself or another image. In the former case, the results are 
exact within numerical noise. In the latter case, because we are correcting an image A with 
a different one B, the method is not exact, but in most cases still provides a quite accurate 
result. 



D. Uncertainty Estimates in Common-Resolution Images 

The x 2 maps also shed light on the estimated flux uncertainties. The value of 



depends not only on the ability of the model to fit the data, but also on the adopted un- 
certainties a\ k . For a given data set and fitted model, underestimation of the uncertainties 
a\ k would lead to higher vales of x 2 ■ Recall that the a\ k were estimated by studying the 
pixel-to-pixel variations in the background regions outside the galaxy mask. If part of the 
background variations actually comes from real variations in emission from foreground dust, 
or dusty background galaxies, then our dust model procedure may actually be able to par- 
tially fit the background variations, leading to a low value of x 2 ■ On the other hand, image 
artifacts and departures of the real data from the models will lead to higher values of % 2 . 
It is clear, for example, that the differences in the PACS and MIPS photometry must lead 
to an increase in the x 2 values, since the model cannot fit simultaneously two different flux 
values. 

If we have N measures from an ideal model with M adjustable parameters, each measure 
having independent gaussian noise with variance a\ , then the x 2 °f the fit will have a chi- 
square distribution with dof degrees of freedom, where dof = N — M. The % 2 /dof map will 
have mean 1, and dispersion ^/2/dof. 

In our background noise estimation procedure, after we recognize and mask all the 
background sources and subtract a "tilted plane" , in each image there remains a zero-mean 
background due to unresolved dim sources. We use the dispersion of the background pixels 
(with respect to the best fit plane) to estimate the stochastic background uncertainty. This 
provides a valid uncertainty estimate since the same background structures will be present 
over the galaxy, so our galaxy fluxes will be uncertain (at least) at this level. However, this 
background dispersion is not only due to random noise in the detectors, it also includes a 
contribution from foreground and background sources. Figure [24] shows the correlation of 




(Dl) 



-79 - 



the dispersion of the background pixels from the MIPS160 (horizontal axis) and SPIRE350 
(vertical axis) cameras. Although the residual distributions have zero mean, there is a strong 
correlation between the MIPS160 and SPIRE350 flux in the background pixels, demostrating 
a component coming from real astrophysical sources. Other cameras show similar correla- 
tions. The modeling can fit these correlated departures from zero by adding positive or 
negative amounts of real dust, reproducing the departures and therefore decreasing the \ 2 
of the fit. We therefore expect to have x 2 /dof < 1 in areas where the flux uncertainties have 
an important contribution from unresolved dim sources (i.e., in the background pixels or low 
surface brightness areas). 



NGC 6946 




'-5 5 10 



MIPS160 (MJy/sr) 

Fig. 24. — Background-subtracted intensities in the MIPS160 and SPIRE350 cameras, for background 
pixels around NGC 6946. There is a strong correlation between the MIPS160 /imand SPIRE350 /im fluxes, 
indicating that the dispersion has a component coming from unresolved dim astrophysical sources. 

Figure 1251 shows the contribution of different cameras to x 2 f° r NCG 628. The fit is done 
at MIPS160 resolution, using all 13 cameras available. We employ the recommended dust 
model with 6 adjusted parameters (fi*, M dust , (?pah, 7, ^min, &) leading to a fit with dof = 7. 
The upper left panel (a) shows the x 2 /dof map. If the data were artifact-free, perfectly 
modeled, with independent gaussian noise, the map would have mean 1, and dispersion 
a/2/7 ~ 0.53. The actual x 2 /dof values are slightly above 1 in the bright regions and below 
1 in the outer parts of the galaxy. In the outer parts of the galaxies, the correlation of 
the noise between the cameras drives the x 2 /dof below the expected value 1. In the bright 
spots of the galaxy, the differences in PACS and MIPS photometry decrease the quality of 
the fit. We estimate the PACS and MIPS photometry uncertainties in a way that their 
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contribution to the x 2 maps is controlled and consistent, but the model will try to fit some 
intermediate (probably non-physical) value which may impair its ability of accurately fitting 
the remaining cameras. 

In order to remedy this situation, uncertainties for the PACS and MIPS photometry 
are estimated by comparing the PACS and MIPS images (see below). The uncertainties are 
slightly overestimated, thus giving more weight to the remaining artifact-free cameras. Figure 
EZST d shows the PACS(S)70/MIPS70 ratio map. Stripes aligned with the MIPS sky scanning 
directions near the galaxy center suggest the presence of MIPS image artifacts. Discrepancies 
in the low surface brightness areas (outer parts of the galaxy) suggest differences in the ability 
of the telescopes to measure the low surface brightness areas correctly. The top right panel 
(c) shows the assigned MIPS70 signal/noise map. The procedure described in Appendix D.l 
assigns large uncertainties to the areas where the PACS and MIPS images have discrepancies, 
leading to reduced S/N. This is a critical step toward getting \ 2 maps that really measure 
the goodness of the model fit, not only the data discrepancies. 

The lower row of Figure m] shows the contribution A% 2 of selected cameras to the total 
X 2 for NGC 628. The lower left panel (d) shows the contribution of MIPS24 camera. We note 
that this contribution is extremely low, under 0.05 for most of the galaxy. This is because the 
dust model can adjust the 24/imflux relatively independently of the other observed bands by 
adjusting 7 and a. The 24/imband is unique in lying a factor of 3 in wavelength away from 
the nearest other bands (8/imand 70/im). Figure [23e is the contribution of MIPS70/imto 
the total x 2 - Even though the MIPS70 camera has severe artifacts, its A% 2 near the center 
is controlled by the large values of a assigned. Figure [25T shows A% 2 for SPIRE350. The 
camera behaves exactly as expected, contributing A% 2 0.6 near the galaxy nucleus and 
smaller A% 2 in the low surface brightness areas. The behavior for NGC 6946 is similar. 

D.l. Uncertainty Estimate Procedure 

In order to estimate uncertainties for the different cameras, we proceed as follows. 

For each camera of wavelength after the image processing (rotation to RA-Dec, 
background subtraction, convolution to a common PSF, correction for boundary effects and 
missing data, and resampling to the final grid) the flux in each final pixel is a (known) linear 
combination of the flux of (in principle) all the original pixels of the camera. If the statistical 
properties of the uncertainties in the original pixel fluxes were known, it would be possible 
to infer the uncertainties (and their statistical properties) in each final pixel. The original 
maps oversample the beam, and have artifacts that extend over several pixels, so realistic 
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Fig. 25. — Analysis of the % 2 contributions for NCG 628. The fit is done at MIPS160 resolution, using 
all 13 cameras with 6 adjusted parameters (fi*, Mdust 5 (7pah, 7, U m - m , a). Top left panel: x 2 /dof map. Top 
center panel: PACS70/MIPS70 ratio map. Top right panel: estimated signal/noise map for MIPS70. The 
lower row shows the contribution Ax 2 of selected cameras to the total % 2; left panel shows MIPS24 camera, 
center panel MIPS70, and right panel shows SPIRE350. The A\ 2 behaves as expected. 
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statistical properties of the uncertainties are difficult to determine. We therefore estimate 
the uncertainties directly in the final (post-processed) image. 

We always assume that the pixel noise is not correlated between different cameras, so 
in what follows we work in each image independently, and omit the subindex identifying the 
camera in the following discussion. For each pixel, we will consider two independent sources 
of uncertainties 5 sto and 5 sys . The term 5 sto will be the uncertainty due to stochastic pixel 
noise and the existence of (unresolved) background dim sources. We estimate it by measuring 
the flux fluctuations in the background pixels. The term 5 sys will include systematic and 
calibration uncertainties. 

In order to estimate the term 5 sto , the stochastic pixel noise, we proceed as follows. We 
re-grid the background mask (BM) to the final-map grid to define a set of iVb g "background 
pixels". For each camera k, we compute the background dispersion S sto ,x k as: 



V { bg ' (x,y)€(BM) 

where 7£ 6s is the background-subtracted flux of the camera of wavelength in the pixel (x, y) 
in the last stage of the image processing (i.e., after the image has been rotated to RA-Dec, 
background plane subtracted, convolved to a common PSF, corrected for boundaries effects 
and missing data, and resampled to the final grid). The background dispersion S sto ,x k includes 
the propagation of original pixel noise and artifacts into the background pixels, possible non- 
ideal background subtraction and the contributions of unidentified faint background sources. 
$sto,\ k does not contain information on the correlation of the noise among the different pixels 
or calibration uncertainties. 

For each camera, we consider the pixel-dependent uncertainty S sys j as 

^=max{6X If 3 ,^} , (D3) 

where e is the calibration uncertainty for the camera; K is for IRAC, MIPS24 and SPIRE 
cameras and an estimate of systematic uncertainties for MIPS70,160 and PACS cameras 
(where such estimation is possible). The main difficulty in estimating e is that the cameras are 
typically calibrated using point sources. We are performing resolved studies of an extended 
sources, so we have additional uncertainties due to the extended source corrections. 

For the IRAC came ras, the standard calibration uncertainty is e = 0.05 for point 



sources (IDale et al.l 120051 ) . but we will adopt a larger value to account for the uncertain- 
ties in the extended emission correction. The IRAC images are calibrated for point sources. 
We multiplied the intensities by the asymptotic (infinite radii) value of the aperture cor- 
rection factors C = 0.91,0.94,0.66,0.74 for the 3.6/im, 4.5/im, 5.8/imand 8.0/mibands 
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respectively. This correction would lead to a correctly calibrated image if the source had 
constant surface brightness. Since galaxies contain nonuniform emission, the aperture cor- 
rection factor should, in principle, be different in each region. We multiply the entire im- 
age by the factors C, and arbitrarily assign 1/3 of the correction as uncertainty. We use 
e = 0.05 + 1/3 x (1/C-l) = 0.083,0.071,0.221,0.167 for the 3.6/xm, 4.5/im, 5.8/imand 
8.0/iinbands, respectively. These corrections will overestimate the uncertainties in the dif- 
fuse regions and underestimate the uncertainties near bright point sources. 

For the MIPS, PACS and the SPIRE cameras we adopt e = 0.10. T his value is larger 
than the estimated 4% ca libration uncertain ty for point sources for MIPS24 ( lEngelbracht et al 
2007h . « 4% for PA CS jMliller et alJboilh . and 7% for SPIRlfl and compara ble to esti- 
mates for MIPS70 JGordon et al.l l2007h and MIPS160 Jstansberrv et al.l l2007h . For the 



MIPS70,160 and PACS cameras, however, 5 sys will be dominated by K in most of the pixels. 

For MIPS70 and MIPS 160, we can estimate the systematic uncertainties K by compar- 
ing the PACS and MIPS resolved photometry. In most of the pixels, the observed differences 
are larger than expected due to the differences in the camera spectral responses or calibration 
uncertainties. For PACS70 and PACS160 we can estimate the systematic uncertainties K 
by comparing the PACS and MIPS resolved photometry For the PACS100 camera we can 
estimate the systematic uncertainties K by extrapolating the uncertainties estimated for the 
PACS70, and PACS160 cameras. We proceed as follows. 



D.l.l. Systematic Uncertainties for MIPS10 and MIPS160 Cameras. 

Since for the 70/imand 160/xm bands the PACS FWHM are smaller than the corre- 
sponding MIPS FWHM (see Tabled]), in all the common-PSF maps where we employ MIPS 
data we have PACS data as well. Therefore, we can compare PACS and MIPS images after 
the image processing. The comparison is straightforward since both images are in the same 
final-map grid and PSF. 

We consider the difference: 

Dpyi,j = I^PACS(S)j — ^MIPSjl > (D4) 

where /mipsj an d ^pacs(s) j are ^ ne observed MIPS and PACS (Scanamorphos) flux, in the 
same band (70/xmor 160/im). We consider the Scanamorphos data reduction more reliable 



See the SPIRE Observers Manual HERSCHEL-DOC-0798, version 2.4 
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than the HIPE, so we only employ the Scanamorphos images in the PACS - MIPS compar- 
ison. The term -DpMj captures both image artifacts and differences in camera performance 
and calibration. Unfortunately, artifacts in the PACS cameras will induce uncertainties in 
the MIPS fluxes and vice-versa, but there is no robust, automatic way of isolating the image 
artifacts. We take 

-Kmipsj = DpM,j ■ (D5) 



D.1.2. Systematic Uncertainties for PACS70 and PACS160 Cameras. 

Based on our dust modeling, we consider the Scanamorphos data reduction to be more 
reliable than the HIPE-only data reduction. When the MIPS and PACS cameras for a band 
(70/umor 160/im) are used in a final-map PSF (e.g., MIPS70 and PACS70 can be used if the 
final- map PSF is SPIRE350, SPIRE500 or MIPS 160) we proceed in a similar way as for the 
MIPS bands, taking 

-Kpacsj = -DpMj • (D6) 

When the MIPS camera for a band (70/imor 160/im ) is not used in a given final-map 
PSF (e.g., MIPS70 cannot be used if the modelling is done at PACS160 PSF), the comparison 
between PACS and MIPS is more complex. We start by convolving the PACS image into 
the corresponding MIPS PSF, and re-binning it into the MIPS native pixel grid, denoted as 
-^pacsom- We can compare /pAcs@M an d ^mTps s i nce both images are expressed at the same 
grid and have the same PSF. We construct epM, the fractional PACS-MIPS difference (with 
respect to the PACS image), at the MIPS PSF, as 

I jobs _ jobs 

l J PACS(S)@M,J J MIPS,J| 

e PM, J ~ VTobS 1 • \P <) 

l J PACS(S)@M,jl 

We re-bin €pm into the final-map grid, and we take 

E PM,j = e PM,J X |/pACS(S)jl > ( D8 ) 

where J is the pixel in the PACS@M image grid that contains pixel j. This procedure 
produces an uncertainty image -Epmj at the final-map grid and PSF, that if degraded into 
the corresponding MIPS PSF, would be ~ I-^pacs(s)@m ~ -^mipsI- I n this situation, we finally 
take 

-Kpacsj = Ep M ,j ■ (D9) 
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D.1.3. Systematic Uncertainties for PACS 100 Camera. 

In the case of PACS100, we do not have any other camera available to perform a 
direct comparison. We therefore estimate the uncertainties by extrapolating the uncertainties 
previously calculated for PACS 70 and PACS 160 as follows. We therefore will assume: 

#PACSiooj = max{i^ 0ji , K' lmj } , (D10) 

where 

tw _ _^pacs7oj | TObs | mm 

""■70 J - TJobs | X l J PACS(S)100,jfl \ ULL ) 

l J PACS(S)70jl 

tw _ -^PAcsieoj | T o b s i mi9^ 

-"-160J — i jobs I l J PACS(S)100,jl ■ l^-^J 

l J PACS(S)160,jl 

Figure [26] shows the resulting S/N estimates for MIPS70 and PACS70 observations of 
NGC 6946. 



E. Uncertainty Estimates for Model Parameters 



In order to estimate the model parameter uncertainties, we generate a set r = 1, 2, N r 
realizations I\ h j >r including random noise, for each pixel j, modeling the noise I\ h) j >r — 
in pixel j as a sum of 3 components: 



jobs 
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\k,j,r,ind 



sys- 
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Fig. 26. — 70Aimimages of NGC 6946 convolved into MIPS160 PSF. Left: map of PACS70 /MIPS70. 
Center: S/N map of MIPS70. Right: S/N map of PACS70. The image artifact that extends from the galaxy 
nucleus toward the WNW induces a low S/N area in the cameras. This artifact, aligned with the MIPS sky 
scanning direction, is probably due to saturation, latency or non-linearities in the MIPS detector. In other 
areas with PACS70/MIPS70 56 1 it is not clear which camera is responsible for the bad behaviour. 
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8I\ k ,j,r,ind represents the non-correlated noise component of the pixel j, while 8I\ h>r>car and 
8I\ k ,j,r,sys represent the correlated noise component and systematic/calibration uncertainties, 
assumed to be (completely) correlated among the different pixels of the camera. SI\ k j, r ,ind is 
drawn from an independent gaussian distribution for each pixel, while SI\ k , r ,cor and SIx kl j >r ,sys 
are drawn from a global gaussian distribution, properly normalized for each pixel. 

For each camera, we generate a set of N x N r "pixel" gaussian variables, with zero mean 
and variance 1, rjj^ r ,j = 1,2, ...N, r = 1,2, ...,N r , where N is the total number of pixels in 
the image, and a set of 2N r "global" gaussian variables, with zero mean and variance 1, 

lfa.,lf>2-,ll>Nr7<fa.><fa, ■••)0JV r - 

We set: 

5h k ,j,r,ind = X 8*o,\ h X Vj,r (E2) 

SI\ k ,r,cor = x $sto,\ k X 1p r (E3) 

3I*k,j,r,sys = $sys,\ k x 4>r , (E4) 

where S sto ,x k and 8 sySj \ k were defined by eq. ID2I and eq. ID3I The division of the dispersion 
$sto,\ k into equal correlated and uncorrelated components SI\ kt j, r ,ind an d 8I\ k , r ,cor is arbitrary, 
but seems to capture the main features of the noise. A more sophisticated treatment would 
use the statistic of the background region to characterize the correlation of the noise in pixels 
as a function of their separation, but this is beyond the scope of the present work. 

For each pixel, the 1-a uncertainty in the measured flux density a\ k j used to calculate 
the Xj °f the model (see equation [13]) is given by: 

(rx k j = ^ (Ssto,x k ) 2 + (6 sys , Xk y . (E5) 

Uncertainties for the inferred (pixel) model parameters are computed using equation [281 In 
order to compute the uncertainties in the global quantities, we use equation [28] on the global 
quantities, i.e., for each random realization we compute the global inferred quantity and 
finally we compute the dispersion on the inferred global quantities. This approach preserves 
the correlation in the pixel noise consistently into the global uncertainties. 

F. Instrument Comparison: PACS-MIPS Photometry disagreement 

In §7 we show that maps made at PACS160 resolution (using IRAC, MIPS24, and PACS 
data only) are less reliable. Furthermore, maps using IRAC and PACS data at MIPS 160 
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resolution, provide less accurate^] results than maps computed with IRAC and MIPS data 
using the same PSF. We therefore proceed to analyze the differences in the PACS and MIPS 
photometry, to shed light onto the roots of such discrepancies. 

Figures [27] and |28] compare the two different PACS data reduction schemes (Scanamor- 
phos and HIPE) with the MIPS camera for NGC 628 and NGC 6946 respectively. All images 
were convolved into a Gaussian PSF with 64" FWHM, a PSF broad enough that any camera 
PSF mismatch should not produce any artifact. The left column show the PACS/MIPS ratio 
maps over the galaxy mask. The center columns show an intensity-intensity scatter plot for 
the galaxy pixels (the pixels in the left column). The right column has a zoom of the -5 
- 30 MJy/sr region of the center column plots. We compare PACS70 Scanamorphos and 
MIPS70, PACS 70 HIPE and MIPS70, PACS 160 Scanamorphos and MIPS160, and PACS160 
HIPE and MIPS160. In the center and right column, the horizontal and vertical red lines 
correspond to the la error estimates. By the inclusion of the PACS(S)-MIPS difference term 
in the uncertainty estimates (see Appendices ID.1.T1 and [TJ. 1.2j) . the la error estimates will 
intersect the PACS=MIPS line. 

The Scanamorphos and HIPE pipelines give similar results in the high surface bright- 
ness areas, having slightly larger intensities than MIPS (see figures |27] and 128]) . The main 
difference between the two data reduction schemes is in the low surface brightness areas: 
Scanamorphos intensities are larger than HIPE intensities. 

Even though the PACS and MIPS cameras have different spectral response, smooth 
spectra appropriate to star-forming galaxies should result in nearly the same measured 
intensity for both MIPS70 and PACS 70, and likewise for MIPS160 and PACS 160. For 
our best-fit SEDs for NGC 628 and NGC 6946, we expect PACS70/MIPS70 -1.085, and 
PACS160/MIPS160 ~ 1.022 (see Table E]). However, the measured PACS global flux densi- 
ties do not agree so well with MIPS values, as seen in Table El The situation is more critical 
in resolved studies. Even though the global PACS and MIPS fluxes agree within « 20%, the 
ratio maps shows strong discrepancies in the three data sets: the band ratios differ from « 1 
in most of the galaxy. 

Figure I2"9l shows the 70/im /160/im color inferred from the different data sets. Again, all 
images were convolved into a Gaussian PSF with 64" FWHM. The top row corresponds to 
NGC 628 and the bottom row to NGC 6946. The left column is the PACS Scanamorphos 
color, PACS HIPE is shown in the center column, and, MIPS is shown in the the right 
column. The 70/im /160/im colors inferred with the different cameras are also quite different, 
especially in the outer parts of the galaxies. 



That is, the resulting dust/H mass ratios seem less "reasonable' 
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Fig. 27. — Comparison of the different 70/imand 160/imdata available for NGC 628, after background 
subtraction, and convolved into a Gaussian PSF with 64" FWHM. The left column shows the PACS - MIPS 
ratio maps over the galaxy mask. The center column shows scatter plot of the camera intensity of the 
galaxy pixels. The right column is a zoom of the -5 - 30 MJy/sr region of the center column plots. Top 
row: PACS 70 Scanamorphos and MIPS70. Second row: PACS70 HIPE and MIPS70. Third row: PACS160 
Scanamorphos and MIPS160. Bottom row: PACS160 HIPE and MIPS160. The ratio maps shows strong 
discrepancies in the three data sets (i.e. the band ratios differ from w 1 in most of the galaxy). 
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Fig. 29. — Comparison of the 70/!m/160/im color inferred from the different data sets. All images were 
previously convolved into a Gaussian PSF with 64" FWHM. The top row corresponds to NGC 628 and the 
bottom row to NGC 6946. Left column: PACS Scanamorphos color. Center column: PACS HIPE. Right 
column: MIPS. The 70/im/160/!m colors inferred with the different cameras are quite different, specially in 
the outer parts of the galaxies. 
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Table 6: PACS70/MIPS70 flux ratio predicted by the model and observed. 



Band ratio 


NGC 628 
Model Observed 


NGC 6946 
Model Observed 


PACS(H)70/MIPS70 
PACS(S)70/MIPS70 


1.083 0.919 
1.083 1.243 


1.088 1.040 
1.088 1.246 


PACS(H)160/MIPS160 
PACS(S)160/MIPS160 


1.024 0.886 
1.024 1.058 


1.020 1.016 
1.020 1.252 



The colors inferred from PACS(S)70/PACS(S)160 and MIPS70/MIPS160 are similar for 
most of the galaxy, and thus the inferred dust parameters (U min ) obtained at PACS160 resolu- 
tion are similar to the ones obtained at MIPS160 resolution. The large PACS(S)160/MIPS160 
flux ratio makes the modeling at PACS160 resolution add more dust over all the galaxy, com- 
pared to our modeling at MIPS160 resolution, when MIPS constraints are present. 

The colors inferred from PACS(H)70/PACS(H)160, and MIPS70/MIPS160 are very 
different in the outer parts of the galaxy, potentially giving different best-fit parameters in 
the dust modeling. The HIPE pipeline uses a high-pass filter in the time-line observing 
series before constructing the two dimensional image map. In studies of HIPE performance 
using simulated data with realistic noise statistics, it was shown that HIPE intensities are 
always equal to or smaller than the "true" intensities in the extended low surface brightness 
regions (Sauvage 2011, private communication). When only one SPIRE band (SPIRE250) is 
included, then the extremely low flux of PACS(H)160 in the low surface brightness regions 
of the galaxy makes the dust appear to be very cold, and thus increases the amount of dust 
needed to reproduce the SPIRE250 flux in those regions. We therefore do not recommend 
to use the HIPE data reduction pipeline for the PACS cameras. 

We do not recommend working at PACS160 resolution, unless there is reason to think 
that high signal/noise has been achieved. We note, however, that even in regions of mod- 
erately high surface brightness, the PACS and MIPS photometry often disagree at the 30% 
level (see Figs. |27] and 128]) . making dust modeling risky if PACS provides the only data 
longward of 24/im. If SPIRE250 data are available, dust modeling can be done at SPIRE250 
resolution with results that are reasonably reliable. We note in Fig. [T4] that total dust 
masses obtained using IRAC, MIPS24, PACS, and SPIRE250 as constraints are in error by 
only PS 35%, provided the PACS-Scanamorphos data are used. 

We further note that NGC 6946 is a particularly challenging galaxy to model: it has 
over 4 decades of dynamic range in the PACS 160 band, including an extremely high surface 
brightness nucleus. 



